GDAL for Python: extracting subdomains from NetCDF file?

Question

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?

• mdsumner: Why reproject? Especially, why reproject to utm? What have you tried already with GDAL? This reads like "here is my task, please do for me"
• Niccolo Bonfadini: Hello and thanks for passing by. I would need UTM because these rasters will be input to an external software that uses UTM. I am new to GDAL. I am wondering if with GDAL it is possible to make things faster than with ArcPy by clipping the netcdf at once at the beginning. I apologize if my question seemed too pretentious. Thanks
• mdsumner: It is actually really straightforward with GDAL. Start by getting the gdalinfo output of your file, and the variable within that you want. Probably it will present as a subdataset in the file.
• nickves: Can you please post your gdalinfo of your netCDF? $ gdalinfo myfile.nc

Answer 1

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.

Answer 2

netCDF4-python will let you subset (using numpy slicing syntax) the data variables without reading the full data from the disk.

• Niccolo Bonfadini: awesome! slicing techniques you mean like this? from numpy import array >>> a = array([[1,2,3],[3,4,5],[4,5,6]]) >>> a[:,1:] array([[2, 3], [4, 5], [5, 6]]) ?? I was also wondering, if I use such technique with numpy arrays then i will lose all the informations contained in the header of the dataset file i get by using the Dataset() method in NETCDF4....???
• DopplerShift: Yes, that slicing will work. When you slice a NetCDF variable, it does become a NumPy array, which does not support all the metadata and attributes that can be attached. However, you can still retrieve this information from the original variable object.