0

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?

4
  • 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"
    – mdsumner
    Commented Sep 3, 2014 at 11:41
  • 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 Commented Sep 3, 2014 at 11:43
  • 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.
    – mdsumner
    Commented Sep 3, 2014 at 11:49
  • Can you please post your gdalinfo of your netCDF? $ gdalinfo myfile.nc
    – nickves
    Commented Sep 3, 2014 at 12:50

2 Answers 2

7

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.

0
0

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

2
  • 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....??? Commented Sep 4, 2014 at 9:34
  • 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. Commented Sep 15, 2014 at 15:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.