Search code examples
pythongisnetcdfqgisnetcdf4

Extract the values of netCDF variables over many time steps at point locations


I have a netCDF file that consists of several meteorological variables (e.g. boundary layer height, surface temperature, surface pressure) measured once daily for a 10-year period over a region. As a result, there are over 4,000 time steps per variable.

I would like to use QGIS or Python to extract the value of each meteorological variable at several point locations for the entirety of the time series. The desired output is either (a) a new vector layer consisting of the point data and the extracted meteorological variable values over all time steps, or (b) a similar output, but in CSV/Excel format.

The "Add raster values to points" tool in QGIS performs the desired operation, but I can't apply this to my NC file since it is mesh data, not raster. Any help here would be appreciated. I'm also quite comfortable with Python is there's a Python package that could achieve this.

Here's the output of ncdump for my netCDF file:

netcdf Beijing_4UTC_2000-10 {
dimensions:
        longitude = 36 ;
        latitude = 19 ;
        time = 4018 ;
variables:
        float longitude(longitude) ;
                longitude:units = "degrees_east" ;
                longitude:long_name = "longitude" ;
        float latitude(latitude) ;
                latitude:units = "degrees_north" ;
                latitude:long_name = "latitude" ;
        int time(time) ;
                time:units = "hours since 1900-01-01 00:00:00.0" ;
                time:long_name = "time" ;
                time:calendar = "gregorian" ;
        short u100(time, latitude, longitude) ;
                u100:scale_factor = 0.000622911066002879 ;
                u100:add_offset = 0.118876376345661 ;
                u100:_FillValue = -32767s ;
                u100:missing_value = -32767s ;
                u100:units = "m s**-1" ;
                u100:long_name = "100 metre U wind component" ;
        short v100(time, latitude, longitude) ;
                v100:scale_factor = 0.000694160650027236 ;
                v100:add_offset = -3.10686248788727 ;
                v100:_FillValue = -32767s ;
                v100:missing_value = -32767s ;
                v100:units = "m s**-1" ;
                v100:long_name = "100 metre V wind component" ;
        short u10(time, latitude, longitude) ;
                u10:scale_factor = 0.000472531631472288 ;
                u10:add_offset = -0.729659788764955 ;
                u10:_FillValue = -32767s ;
                u10:missing_value = -32767s ;
                u10:units = "m s**-1" ;
                u10:long_name = "10 metre U wind component" ;
        short v10(time, latitude, longitude) ;
                v10:scale_factor = 0.000513328716515632 ;
                v10:add_offset = -3.10863126000035 ;
                v10:_FillValue = -32767s ;
                v10:missing_value = -32767s ;
                v10:units = "m s**-1" ;
                v10:long_name = "10 metre V wind component" ;
        short d2m(time, latitude, longitude) ;
                d2m:scale_factor = 0.00102728301940017 ;
                d2m:add_offset = 268.55166378769 ;
                d2m:_FillValue = -32767s ;
                d2m:missing_value = -32767s ;
                d2m:units = "K" ;
                d2m:long_name = "2 metre dewpoint temperature" ;
        short t2m(time, latitude, longitude) ;
                t2m:scale_factor = 0.00110899471613419 ;
                t2m:add_offset = 279.266550605181 ;
                t2m:_FillValue = -32767s ;
                t2m:missing_value = -32767s ;
                t2m:units = "K" ;
                t2m:long_name = "2 metre temperature" ;
        short blh(time, latitude, longitude) ;
                blh:scale_factor = 0.0884957109728695 ;
                blh:add_offset = 2911.24127842015 ;
                blh:_FillValue = -32767s ;
                blh:missing_value = -32767s ;
                blh:units = "m" ;
                blh:long_name = "Boundary layer height" ;
        short sp(time, latitude, longitude) ;
                sp:scale_factor = 0.374649985503487 ;
                sp:add_offset = 92661.2189250073 ;
                sp:_FillValue = -32767s ;
                sp:missing_value = -32767s ;
                sp:units = "Pa" ;
                sp:long_name = "Surface pressure" ;
                sp:standard_name = "surface_air_pressure" ;
        short tp(time, latitude, longitude) ;
                tp:scale_factor = 3.21524733289788e-07 ;
                tp:add_offset = 0.0105350794109732 ;
                tp:_FillValue = -32767s ;
                tp:missing_value = -32767s ;
                tp:units = "m" ;
                tp:long_name = "Total precipitation" ;

// global attributes:
                :Conventions = "CF-1.6" ;
                :history = "2020-06-30 13:19:05 GMT by grib_to_netcdf-2.16.0: /opt/ecmwf/eccodes/bin/grib_to_netcdf -S param -o /cache/data3/adaptor.mars.internal-1593521146.6122832-30945-36-b3534028-e955-439d-b3fd-a960d29aa305.nc /cache/tmp/b3534028-e955-439d-b3fd-a960d29aa305-adaptor.mars.internal-1593521146.6128585-30945-7-tmp.grib" ;
}

Solution

  • Managed to solve this using the 'Make NetCDF Table View" tool in ArcGIS Pro. I selected all the variables of interest, set the row dimension to the 'Time' variable, and used lat/long coordinates as dimension-value pairs. This returns a table with the value of the variable (e.g. surface temperature) for each time step in the netCDF file sampled at the specified XY location.