Zonal statistics with GDAL only (not Python)

Question

I'm helping an undergrad with her dissertation. She has monthly mean data from ERA-5 for three variables in netcdf format for the past 20 years. She wants a CSV or Excel file which contains the zonal statistics for the shapefile covering her region only for the three variables for four statistical tests: mean, min, max and variance. Each row would have a date column, then 12 further columns with mean/min/max/variance for the three variables. Thus, 13 columns. Normally I capture the GDAL command from QGIS and compose a Linux shell script (I am very experienced writing bash shell scripts including awk and grep - I don't need help with this).

I don't see a non Python-GDAL solution for calculating these quantities. Ideally I'd like to avoid using Python code in favour of a GDAL command-line-only command I could use in a Linux shell script -- a GDAL command which takes as input the netcdf file (the raster in this instance) and the region of interest as a shapefile, then outputs four statistical quantities for me to capture to place into a CSV file using awk/grep/etc.

Does this functionality exist in non-Python GDAL?

Answer 1

There are a bunch of open source tools beyond the GDAL NetCDF driver you could use. These tools are specialized for your outlined problem area of weather/climate statistics.

1. You could have a look "Climate Data Operators" to work with NetCDf /GRIB encoded data of the ERA-5 ECMWF reanalysis data. CDO is focused on "climate data" with as well defined interface. This interface addresses an international standard with a specialized NetCDF/GRIB data configuration for provenience, dimensions (spatial, time) and measured/calculated records (pressure, temperature ...). AFAIK most ECMWF data products follow this standard due to the need of international exchange of data sets in the context of (European Center For) Medium-Range Weather Forcast. The toolkit is very fast and has an operator set for statistics -- section 1.6.4 -- page 89 ff.. For every latitude is an operator over all longitudes defined. Operators for zonal statistics are:

       cdo [zonemin,zonmax,zonmean,zonmean,zonvar,zonstd] input-file output-file
       cdo zonepctl,num input-file output-file

2. Another toolkit is "NetCDF operators". To get data out of the NetCDF result set, you could have a look at "How can I extract data from a NetCDF file for a specific location?". NCO has also a set of statistic operators -- section 4.6 -- page 254 ff.. The syntax for an aggregation (min, max, avg, stddev, etc.) of data records over an dimension (spatial i.e. longitude, latitude and time...) is more complex, because you can use a free NetCDF format to store the matter.

    nces [-3] [-4] [-5] [-6] [-7] [-A] [-C] [-c] [--cb y1,y2,m1,m2,tpd]
    [--cnk_byt sz_byt] [--cnk_csh sz_byt] [--cnk_dmn nm,sz_lmn]
    [--cnk_map map] [--cnk_min sz_byt] [--cnk_plc plc] [--cnk_scl sz_lmn]
    [-D dbg] [-d dim,[min][,[max][,[stride]]] [-F]
    [-G gpe_dsc] [-g grp[,...]] [--glb ...]
    [-h] [--hdf] [--hdr_pad nbr] [--hpss]
    [-L dfl_lvl] [-l path] [-n loop]
    [--no_cll_msr] [--no_frm_trm] [--no_tmp_fl] [--nsm_fl|grp] [--nsm_sfx sfx]
    [-O] [-o output-
    file] [-p path] [--ppc ...] [-R] [-r] [--ram_all] [--rth_dbl|flt]
    [-t thr_nbr] [--unn] [-v var[,...]] [-w wgt] [-X ...] [-x] [-y op_typ]
    [input-files] [output-file]

3. For completeness I will mention the NetCDF Command Language. You could write scripts in NCL and have predefined functions for statistics too. NCL produces terrific predefined maps. These maps are good, if you need an orientation and the use cases are well documented. The example page for zonal statistics provides eight use cases. NCL example for a zonal average temperature (www.ncl.ucar.edu).

    ;*************************************************
    ; zonal_3.ncl
    ;
    ; Concepts illustrated:
    ;   - Attaching a zonal means plot to a cylindrical equidistant contour plot
    ;************************************************
    ;
    ; These files are loaded by default in NCL V6.2.0 and newer
    ; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
    ; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
    ;************************************************
    begin
    ;************************************************
    ; variable and file handling
    ;************************************************
      in  = addfile("83.nc","r") 
      ts = in->TS                                   ; select variable to ave
    ;************************************************
    ; plotting
    ;************************************************
      wks  = gsn_open_wks("png","zonal")            ; send graphics to PNG file
      
      res                 = True                ; make plot mods
      res@tiMainString    = "Zonal Average"         ; Title for the plot
      res@gsnZonalMean    = True                    ; put zonal on side
      res@gsnZonalMeanXMinF = 240.          ; set minimum X-axis value for zonal mean plot  
      res@gsnZonalMeanXMaxF = 315.          ; set maximum X-axis value for zonal mean plot  
      res@gsnZonalMeanYRefLine = 273.15     ; set reference line X-axis value
      
      plot=gsn_csm_contour_map(wks,ts(0,:,:),res)   ; plot temp contours 
    ;***********************************************
    end

Answer 2

You can use pkextract. Requires an image and vector layer (point, line, polygon), and you can decide whether to calculate a value from the pixels or extract a random point.

Usage: pkextract -i input.tif -s sample.shp -o output.shp -r median

Details: http://www.nongnu.org/pktools/html/md_pkextract.html

Examples: http://www.nongnu.org/pktools/html/md_examples_pkextract.html#examples_pkextract

Install: http://www.nongnu.org/pktools/html/index.html#pktools_installation

Answer 3

This is not a command-line solution however given the nature of need in the question, I think it might help be helpful.

If my understanding is correct, you want to get the summary statistics for the shapes of a given shapefile from a netcdf file (for example temperate for sub-basins of a larger basin). This will include two steps of remapping the variables to the shapefile and then extracting the summary statistics such as min or max.

I would suggest using the easymore python package. What you need to do is just to replace the netcdf files, shapefile, and the variables based on the link to the example for era5. The remapped output is areal average of the variables for each shape. The final netcdf remapped file(s) is(are) a netcdf with dimension of time (at original/source netcdf) and number of shapes in the shapefile. After this step, the statistical analysis can be done on the remapped variables for min, mean, or max or any other summary statistics.