Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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?

Any feedback will be greatly appreciated!!

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

Edit

Hi Luke (Not sure how to post this as a comment.)

I have been able to convert my shapefile to a raster using gdal_RasterizeLayer link you 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. Thanks.

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
share|improve this question
    
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. –  Luke Nov 16 '11 at 1:33
    
@questions I've updated my answer to include an example using your extents. –  Luke Nov 17 '11 at 1:41

2 Answers 2

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.

share|improve this answer

You are looking for the gdal.RasterizeLayer function.

Below, if it is of any use, is a script I used to use before the gdal_rasterize utility was able to create the output raster on the fly. Apologies for the mass of code but it does show you various ways of using gdal.RasterizeLayer.

'''
Usage:
       gdal_rasterize.py [-p prototype_raster] [-u]
       [-at] [-burn value] | [-a attribute]
       [-where expression] | [-sql select_statement]
       [-of format] [-a_srs srs_def] [-co \"NAME=VALUE\"]*
       [-a_nodata value] [-init value]
       [-te xmin ymin xmax ymax] [-tr xres yres] [-ts width height]
       [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/
            CInt16/CInt32/CFloat32/CFloat64}]
       <src_datasource> <dst_filename>
Where:
    -p          Path to a prototype raster
    -u          Update the existing destination raster
    -at         Enables the ALL_TOUCHED rasterization option so that
                all pixels touched by lines or polygons will be
                updated not just those one the line render path, or
                whose center point is within the polygon. Defaults to
                disabled for normal rendering rules.
    -burn       A fixed value to burn into a band for all objects.
    -a          Identifies an attribute field on the features to be
                used for a burn in value.
    -where      An optional SQL WHERE style query expression to be
                applied to select features to burn in from the input
                layer.
    -sql        An SQL statement to be evaluated against the
                datasource to produce a virtual layer of features to
                be burned in.
    -of         The output format. Default is GeoTIFF (GTIFF). Use the
                short format name.
    -a_nodata   Assign a specified nodata value to output bands.
    -init       Pre-initialize the output image with this value.
                However, it is not marked as the nodata value in the output
                file.
    -a_srs      Override the projection for the output file. If not specified,
                the projection of the input vector file will be used if
                available. If incompatible projections between input and output
                files, no attempt will be made to reproject features. The
                srs_def may be any of the usual GDAL/OGR forms, complete WKT,
                PROJ.4, EPSG:n or a file containing the WKT.
    -co         Passes a creation option to the output format driver. Multiple
                -co options may be listed. See format specific documentation
                for legal creation options for each format.
    -te         set georeferenced extents. The values must be expressed in
                georeferenced units. If not specified, the extent of the output
                file will be the extent of the vector layers.
    -tr         set target resolution. The values must be expressed in
                georeferenced units. Both must be positive values.
    -tap        Align the coordinates of the extent of the output file to the
                values of the -tr, such that the aligned extent includes the
                minimum extent.
    -ts         Set output file size in pixels and lines. Note that -ts cannot
                be used with -tr
    -ot         For the output bands to be of the indicated data type.
                Defaults to Float64
    -w          Intermediate working format for rasters.
                Defaults to 'MEM'.
'''
#---Imports
from osgeo import gdal, gdalconst, ogr
import sys, os,  math, tempfile

gdal.UseExceptions()
ogr.UseExceptions()

def create_raster(xsize, ysize, driver='MEM',tmpfile='', gt=None, srs_wkt=None, nodata=None,  init=None, datatype=gdal.GDT_Byte):
    # Create a memory raster to rasterize into.
    out_ds = gdal.GetDriverByName(driver).Create(tmpfile, xsize, ysize, 1 ,datatype)
    if init is not None:out_ds.GetRasterBand(1).Fill(init)
    if nodata is not None:out_ds.GetRasterBand(1).SetNoDataValue(nodata)
    if gt:out_ds.SetGeoTransform(gt)
    if srs_wkt:out_ds.SetProjection(srs_wkt)
    return out_ds

def rasterize(src,dst,update=False,
              prototype=None,
              nodata=None,init=None,
              te=None,tr=None, ts=None,
              sql=None, where=None,
              out_format='GTIFF', out_type='Byte',
              out_srs=None, co=[], working_format='MEM',
              **kwargs):

    #Some examples
    #dstds=rasterize(src ,dst,ts=[500, 500], out_type='Int16',  burn_values=[37])
    #dstds=rasterize(src ,'',prototype=prot, nodata=-999, init=-999, ts=[500, 500], out_type='Int16', options=["ATTRIBUTE=test"],  out_format='MEM')
    #rasterize(src ,dst, tr=[0.01, 0.01], nodata=-999, init=-999, out_type='Int16', options=["ATTRIBUTE=test"])
    #rasterize(src ,dst, tr=[0.01, 0.01], options=["ATTRIBUTE=test"],  out_format='ECW',  co=['TARGET=10'])
    #rasterize(src ,dst, where='test=1', update=True, burn_values = [10])
    #rasterize(src ,dst, tr=[0.01, 0.01], te=[143.5, -43.5, 148.5, -39.5], out_type='Int16', options=["ATTRIBUTE=test"])

    #Temporary files
    if working_format=='MEM':tmpfd,tmpfile=[None,'']
    else:tmpfd,tmpfile=tempfile.mkstemp()

    #Open the vector layer to rasterize from
    src_ds = ogr.Open(src)
    if not src_ds: raise RuntimeError,'\'%s\' does not exist in the file system.' % src
    if sql:src_lyr=ds.ExecuteSQL(sql)
    else:
        src_lyr=src_ds.GetLayer()
        if where:src_lyr.SetAttributeFilter(where)
    xmin,xmax,ymin,ymax=src_lyr.GetExtent()
    src_ext=xmin,ymin,xmax,ymax
    try:
        src_wkt=src_lyr.GetSpatialRef().ExportToWkt()
        if not out_srs:out_srs=src_wkt
    except:
        out_srs='LOCAL_CS["arbitrary"]' # From http://svn.osgeo.org/gdal/trunk/autotest/alg/rasterize.py
    datatype=gdal.GetDataTypeByName(out_type)

    #Get a GDAL Dataset to rasterize into
    dst_driver=None
    if update: #Can we update an existing raster
        dstds = gdal.Open(dst, gdalconst.GA_Update)
    else:
        gdal.ErrorReset()
        if prototype:
            protds=gdal.Open(prototype)
            dst_driver=protds.GetDriver()
            dstds = create_raster(protds.RasterXSize, protds.RasterYSize, working_format,
                                  tmpfile,protds.GetGeoTransform(), protds.GetProjection(),
                                  datatype=datatype, nodata=nodata, init=init)
            del protds
        elif te and tr:
            out_gt=(te[0],tr[0],0,te[3],0, -tr[1])
            out_cols=int(math.ceil((te[2]-te[0])/tr[0]))
            out_rows=int(math.ceil((te[3]-te[1])/tr[1]))
            dstds = create_raster(out_cols, out_rows, working_format, tmpfile,
                                  out_gt, out_srs, datatype=datatype, nodata=nodata,
                                  init=init)
        elif te and ts:
            xres=(te[2]-te[0])/float(ts[0])
            yres=(te[3]-te[1])/float(ts[1])
            out_gt=(te[0],xres,0,te[3],0, -yres)
            dstds = create_raster(ts[0], ts[1], working_format, tmpfile,
                                  out_gt, out_srs, datatype=datatype, nodata=nodata,
                                  init=init)
        elif ts:
            xres=(src_ext[2]-src_ext[0])/float(ts[0])
            yres=(src_ext[3]-src_ext[1])/float(ts[1])
            out_gt=(src_ext[0],xres,0,src_ext[3],0, -yres)
            dstds = create_raster(ts[0], ts[1], working_format, tmpfile,
                                  out_gt, out_srs, datatype=datatype, nodata=nodata,
                                  init=init)
        elif tr:
            out_cols=int(math.ceil((src_ext[2]-src_ext[0])/float(tr[0])))
            out_rows=int(math.ceil((src_ext[3]-src_ext[1])/float(tr[1])))
            out_gt=(src_ext[0],tr[0],0,src_ext[3],0,-tr[1])
            dstds = create_raster(out_cols, out_rows, working_format, tmpfile,
                                  out_gt, out_srs, datatype=datatype, nodata=nodata,
                                  init=init)
        else:
            raise RuntimeError, '%s does not exist and neither a prototype raster nor appropriate options were specified!'%dst

        if dst_driver  is None:
            try:dst_driver = gdal.GetDriverByName(out_format)
            except:raise RuntimeError,'Format driver %s not found, pick a supported driver.' % out_format

    err = gdal.RasterizeLayer(dstds, [1], src_lyr, **kwargs)
    if err != 0:
        raise RuntimeError,  "error rasterizing layer: %s" % err
    if not update:
        dstds=dst_driver.CreateCopy(dst,dstds, 1, co)
        try:
            os.close(tmpfd)
            os.unlink(tmpfile)
        except:pass

    dstds.FlushCache()
    return dstds

if __name__ == '__main__':
    try:
        kwargs={}
        args=map(str,  sys.argv[1:])
        src,dst=args[-2:];del args[-2:]
        while args: #This is used instead of getopt or optparse
                    #to support gdal non standard long option name
                    #format "-option value" instead of "--option=value"
            arg=args.pop(0)
            if arg == '-p':
                kwargs['prototype']=args.pop(0)
            elif arg == '-u':
                kwargs['update']=True
            elif arg == '-sql':
                kwargs['sql']=args.pop(0)
            elif arg == '-w':
                kwargs['where']=args.pop(0)
            elif arg == '-at':
                kwargs['options'].append('ALL_TOUCHED=TRUE')
            elif arg == '-burn':
                kwargs['burn_values']=[int(args.pop(0))]
            elif arg == '-a':
                if kwargs.has_key(options):kwargs['options'].append()
                else:kwargs['options']=['ATTRIBUTE=%s'%args.pop(0)]
            elif arg == '-a_srs':
                kwargs['out_srs']=args.pop(0)
            elif arg == '-co':
                if kwargs.has_key(co):kwargs['co'].append(args.pop(0))
                else:kwargs['co']=[args.pop(0)]
            elif arg == '-a_nodata':
                kwargs['nodata']=args.pop(0)
            elif arg == '-init':
                kwargs['init']=args.pop(0)
            elif arg == '-te':
                kwargs['te']=map(float, args[:4]);del args[:4]
            elif arg == '-tr':
                kwargs['tr']=map(int, args[:2]);del args[:2]
            elif arg == '-ts':
                kwargs['ts']=map(int, args[:2]);del args[:2]
            elif arg == '-of':
                kwargs['out_format']=args.pop(0)
            elif arg == '-ot':
                kwargs['out_type']=args.pop(0)
            elif arg == '-w':
                kwargs['working_format']=args.pop(0)
            else:raise Exception
    except:
        print >> sys.stderr,__doc__
        sys.exit(1)

    dstds=rasterize(str(src) ,str(dst),  **kwargs)

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

#If gdal can read the netcdf file directly
mask_ds=rasterize(shapefile ,'',prototype='E:/GriddedData/T_SFC/IDZ00026_VIC_ADFD_T_SFC.nc',  out_type='Byte',init=0, burn_values=[1],  out_format='MEM')
#Else
mask_ds=rasterize(shapefile ,'',te=[xmin, ymin, xmax, ymax], tr=[xres, yres], out_type='Byte',init=0, burn_values=[1],  out_format='MEM')
mask_arr=mask_ds.GetRasterBand(1).ReadAsArray()

EDIT: 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()
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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