*! version 3.0 10oct2025 (unified GeoTIFF/NetCDF)

Title

zonalstats -- Zonal statistics command for GeoTIFF and NetCDF raster files

Syntax

For GeoTiff files:

zonalstats rasterfilename using shapefile, stats(string) [band(#) clear crs(string)]

For NetCDF files:

zonalstats rasterfilename using shapefile, stats(string) var(string) [clear origin(numlist) size(numlist) crs(string)]

rasterfilename can be a GeoTIFF (.tif or .tiff) or NetCDF (.nc) file. The command automatically detects the file type and uses the appropriate processing method.

shapefile must include accompanying .shx .dbf (and ideally .prj) files.

Description

zonalstats computes statistics of raster cell values aggregated over polygon zones defined in a shapefile. It supports both GeoTIFF and NetCDF raster files:

• GeoTIFF files: Direct reading of one band from GeoTIFF files on disk.

• NetCDF files: Reading and processing of NetCDF variables with optional slicing.

The command automatically detects the file type based on the file extension and uses the appropriate processing method.

Features:

• Reprojects the shapefile to match the raster CRS if needed.

• Supports multiple statistics (count avg min max std sum).

• Uses GeoTools Java libraries for spatial processing.

• Automatic CRS detection with fallback to user-specified CRS.

• NetCDF slicing support with origin() and size() options.

Dependencies

The crsconvert command requires Java libraries from GeoTools and netCDF-Java.

Run geotools_init to configure the GeoTools library path.

Run netcdf_init to configure the netCDF-Java library path (pointing to netcdfAll-5.9.1.jar).

Options

Common options

stats(string) Statistics to compute; default avg. Any space separated subset of count avg min max std sum. Invalid names produce an error.

clear Clear current data in memory before loading results (required if data present).

crs(string) Coordinate reference system for the raster data. If the raster file contains CRS information, this option is ignored and a message is displayed.

GeoTIFF-specific options

band(#) Band index (1-based) for multi-band GeoTIFF. Default 1; must be >=1.

NetCDF-specific options

var(string) Variable name in the NetCDF file to process (required for NetCDF files).

origin(numlist) Origin coordinates (1-based) for slicing the NetCDF variable. Must be integers >0.

size(numlist) Size of each dimension for slicing. At most 2 dimensions can have size >1 (2D grid requirement).

Remarks

File type detection. The command automatically detects whether the input file is a GeoTIFF (.tif/.tiff) or NetCDF (.nc) file based on the file extension and uses the appropriate processing method.

NetCDF processing. For NetCDF files, the command supports multi-dimensional data with automatic detection of 2D spatial grids (allowing singleton dimensions). Use origin() and size() options for slicing large datasets. At most 2 dimensions can have size >1 to maintain 2D grid requirements.

CRS handling. The command attempts to automatically detect CRS from raster files. If CRS is found, user-provided crs() is ignored with a notification. If no CRS is detected and crs() is not provided, an error occurs. Shapefile is reprojected to match raster CRS when needed.

Performance tips.

• Limit requested statistics to those needed.

• Use appropriate CRS in projected meters for large area analyses.

• Pre-filter polygons to study region before running.

• For large NetCDF files, use slicing with origin() and size() to process subsets.

• For very large rasters consider tiling externally; current command reads full needed extent.

Examples

Nighttime lights statistics by city (sum + average):

. zonalstats DMSP-like2020.tif using hunan.shp, stats("sum avg") clear

NetCDF with slicing (subset of data):

. local url = "https://nex-gddp-cmip6.s3-us-west-2.amazonaws.com/NEX-GDDP-CMIP6/BCC-CSM2-MR/ssp245/r1i1p1f1/tas/tas_day_BCC-CSM2-MR_ssp245_r1i1p1f1_gn_2050.nc"

. ncread lon using `url'

. gen n=_n

. qui sum n if lon>=108 & lon<=115

. local lon_start = r(min)

. local lon_count = r(N)

. ncread lat using `url', clear

. gen n=_n

. qui sum n if lat>=24 & lat<=31

. local lat_start = r(min)

. local lat_count = r(N)

. zonalstats `url' using "hunan.shp", var(tas) stats(avg) origin(1 `lat_start' `lon_start') size(1 `lat_count' `lon_count') crs(EPSG:4326) clear

Stored results

Output dataset (or results frame in vector mode) contains:

Variable Description

zone_id_vars All non-geometry attributes from shapefile polygons

count Pixel count in zone (if requested)

avg Mean cell value (if requested)

min Minimum cell value (if requested)

max Maximum cell value (if requested)

std Standard deviation (if requested)

sum Sum of cell values (if requested)

Each observation = one polygon zone. Statistics will not include nodata pixels.

Author

Kerry Du

School of Management, Xiamen University, China

Email: kerrydu@xmu.edu.cn

Chunxia Chen

School of Management, Xiamen University, China

Email: 35720241151353@stu.xmu.edu.cn

Yang Song

School of Economics, Hefei University of Technology, China

Email: ss0706082021@163.com

Ruipeng Tan

School of Economics, Hefei University of Technology, China

Email: tanruipeng@hfut.edu.cn

Also see

Online: geotools_init, netcdf_init