Python GDAL – How to Extract Subdomains from NetCDF File Using GDAL for Python

batchcoordinate systemgdalnetcdfpython

I need to make a script that processes a NetCDF file that contains 3 days of hourly forecast data from the norwegian meteorological office.

  1. The NetCDF file contains various data I need (Precipitation,Tmperature,Wind etc).

  2. The NetCDF file is in a lambert projection while I will need to project it into WGS84 UTM 32N.

  3. Also I will need to resample from 2.5km (forecast inputs) to 1km(output) grid cells.

  4. I need to save it into the IDRISI format .rst

PROBLEM! The original NetCDF is HUGE, covering the whole of scandinavia + neighbouring countries. Thus I will need a system that processes quick.

I managed to do this already with ArcPy, but the process was too slow since for every hour timestep I needed to extract one by one the huge rasters, and only then could I clip them down.

Maybe in GDAL (in Python) there is a way to first clip at once the whole netcdf and then continue the processing with a smaller netcdf?

Best Answer

The Norwegian Met office has a THREDDS server at http://thredds.met.no/thredds/ so if you see the forecast you are trying to access there, you can extract just the subset you want from the OPeNDAP URL, which NetCDF4-Python treats like a local netcdf file.

For example:

import netCDF4
url = 'http://thredds.met.no/thredds/dodsC/arome25/arome_norway_default2_5km_latest.nc'

nc = netCDF4.Dataset(url)
ncv = nc.variables
ncv.keys()

# subset and subsample
lon = ncv['longitude'][10:-10:2,20:-10:2]
lat = ncv['latitude'][10:-10:2,20:-10:2]
# read the 1st time step
itime = 0
tair = ncv['air_temperature_2m'][itime,10:-10:2,20:-10:2]

pcolormesh(lon,lat,tair);
colorbar();

produces this plot:

enter image description here

You could process the subset/subsample the data this way, or you could also use NCO tools to subset the OPeNDAP url:

ncks -d x,10,30 -d y,20,40 -d time,0 http://thredds.met.no/thredds/dodsC/arome25/arome_norway_default2_5km_latest.nc  subset.nc

From there you should either be able to use GDAL to convert or use the pyproj with the proj4 parameters included in the file to convert to whatever you need.

Related Question