7

I have been using Python to process huge arrays stored as NetCDF files. I would like to calculate the mean of an area defined by a shapefile. I have just installed GDAL but if there are other tools I should use please let me know.

So my main question is how can I turn my shapefile into a mask?

Then use my mask to calculate the NumPy mean of my array within the area of the mask?

shp="country.shp"
array=mynumpyarray
mask=useGDALtochangeshapefiletomask
meanofarea=N.mean(array,limits=mask)

Edit

I have been able to convert my shapefile to a raster using gdal_RasterizeLayer link supplied and then a masked array using

maskarray=mask_ds.GetRasterBand(1).ReadAsArray()

Now it is just a question of matching my NetCDF extent with my raster/maskarray. Here are my extents as requested.

my shapefile/raster extent: x_min, x_max, y_min, y_max 140.962408758 149.974994992 -39.1366533667 -33.9813898583

my netcdf file extent: min longitude: 139.8 max longitude: 150.0 min latitude -39.2 max latitude: -33.6 LAT size 106 LON size 193
  • Thanks for the feedback but I can't seem to get it to work, mostly because I do not understand it. I have tried running the code from the link you posted for gdal.RasterizeLayer but keep getting an error Traceback (most recent call last): File "d:/documents and settings/gdal_rasterizelayer2.py", line 46, in <module> rasterize(pixel_size=25) File "d:/gdal_rasterizelayer2.py", line 13, in rasterize source_layer = source_ds.GetLayer(0) AttributeError: 'NoneType' object has no attribute 'GetLayer' – questions Nov 16 '11 at 0:48
  • If you can tell me what the extent (xmin,ymin,xmax,ymax) and cellsize of the netcdf raster is (as a comment, not an answer) I will write some code that should just work for you. I do recommend that you look over the code I posted in my answer and try to understand what it's doing though. Basically you need to pass an existing open dataset handle to gdal.RasterizeLayer, most of the code up to that point is designed to create the destination dataset from scratch. – user2856 Nov 16 '11 at 1:33
  • @questions I've updated my answer to include an example using your extents. – user2856 Nov 17 '11 at 1:41
8

You are looking for the gdal.RasterizeLayer function.

You could then use ReadAsArray to turn the rasterized polygon into a numpy array.

Based on your NetCDF file extent and rows/columns, the following code should generate you a numpy 0-1 mask that matches the NetCDF exactly.

shapefile=r'whatever your shapefile path is'
xmin,ymin,xmax,ymax=[139.8,-39.2,150.0,-33.6] #Your extents as given above
ncols,nrows=[193,106] #Your rows/cols as given above
maskvalue = 1

xres=(xmax-xmin)/float(ncols)
yres=(ymax-ymin)/float(nrows)
geotransform=(xmin,xres,0,ymax,0, -yres)

src_ds = ogr.Open(shapefile)
src_lyr=src_ds.GetLayer()

dst_ds = gdal.GetDriverByName('MEM').Create('', ncols, nrows, 1 ,gdal.GDT_Byte)
dst_rb = dst_ds.GetRasterBand(1)
dst_rb.Fill(0) #initialise raster with zeros
dst_rb.SetNoDataValue(0)
dst_ds.SetGeoTransform(geotransform)

err = gdal.RasterizeLayer(dst_ds, [maskvalue], src_lyr)
dst_ds.FlushCache()

mask_arr=dst_ds.GetRasterBand(1).ReadAsArray()
4

Perhaps the easiest solution is to call gdal_rasterize from within your Python code (use subprocess.call or subprocess.check_call), either using your NetCDF file as the destination, or create a separate image file (GeoTIFF is always a good bet) and load it into a numpy array. It may be possible to use GDAL's In Memory Raster format, but I'm not sure what support there is for it in Python, and how it would play with numpy arrays.

There is an overhead with calling a subprocess and loading an image, but that will probably be amortized by the time it takes to do the actual rasterization.

If you want to roll your own, you could write a simple scanline renderer, but I suspect that would be more fraught and possibly slower if written in pure Python.

2

Now in 2017, I use rasterio.mask.mask and fiona to rasterize the shape file and then apply that to the netCDF dataset using xarray.

EDIT: Now in 2019, there is https://github.com/corteva/rioxarray to ease the integration of rasterio and xarray. Clipping a dataset (e.g. from NetCDF) with a geometry is listed as example here: https://corteva.github.io/rioxarray/html/examples/clip_geom.html

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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