0

Problem

I've recently noticed that the output from the gdal.RasterizeLayer function differs from the result I'm getting using the gdal_rasterize utility.

I have geotiff of a SAR product that has been resampled to 500m and my aim is to create a corresponding land mask geotiff. The extent of the data is centered around Florence, Italy and covers the width of Italy with some sea to the east and west.

Sentinel 1A IW VV - Florence, Italy

I'm using the GSHHS L1 vector dataset for the masking and I've written a Python function to perform the operation:

import os
import numpy as np
from osgeo import gdal, ogr, osr

# Source raster details
# src_rast_gt = (597878.6395908913, 500.0, 0.0, 4952884.831090534, 0.0, -500.0)
# src_rast_xsize = 560
# src_rast_ysize = 364
# src_dtype = gdal.GDT_Float32
# src_epsg = 32632

def create_land_mask(src_rast,land_vector='GSHHS_f_L1.shp'):
    # Get raster metadata
    rast_ds = gdal.Open(src_rast)
    if rast_ds is None:
        sys.exit('Failed to open src_file "{}"'.format(src_file))
    band = rast_ds.GetRasterBand(1)
    geotrans = rast_ds.GetGeoTransform()
    xsize, ysize = band.XSize, band.YSize

    mask_ds = ogr.Open(land_vector)
    if rast_ds is None:
        sys.exit('Failed to open land_vector "{}"'.format(land_vector))
    mask_lyr = mask_ds.GetLayer()

    # Create SpatialReference objects to hold raster and vector srs
    rast_srs = osr.SpatialReference(wkt=rast_ds.GetProjection())
    mask_srs = mask_lyr.GetSpatialRef()

    # Get raster corner points
    rast_tl = (geotrans[0] + 0*geotrans[1] + 0*geotrans[2],geotrans[3] + 0*geotrans[4] + 0*geotrans[5])
    rast_tr = (geotrans[0] + xsize*geotrans[1] + 0*geotrans[2],geotrans[3] + xsize*geotrans[4] + 0*geotrans[5])
    rast_br = (geotrans[0] + xsize*geotrans[1] + ysize*geotrans[2],geotrans[3] + xsize*geotrans[4] + ysize*geotrans[5])
    rast_bl = (geotrans[0] + 0*geotrans[1] + ysize*geotrans[2],geotrans[3] + 0*geotrans[4] + ysize*geotrans[5])
    coords = [rast_tl,rast_tr,rast_br,rast_bl,rast_tl]

    # Generate raster extent geometry
    rast_ext = ogr.Geometry(ogr.wkbPolygon)
    lin_ring = ogr.Geometry(ogr.wkbLinearRing)
    for c in coords:
        lin_ring.AddPoint_2D(*c)
    rast_ext.AddGeometry(lin_ring)
    rast_ext.AssignSpatialReference(rast_srs)

    if not rast_srs.IsSame(mask_srs):
        rast_ext.TransformTo(mask_srs)

    # Apply a spatial filter to the shapefile layer.
    mask_lyr.SetSpatialFilter(rast_ext)

    # Create raster to store mask
    drv = gdal.GetDriverByName('GTiff')
    out_path = os.path.splitext(src_rast)[0]+'_mask.tif'
    if os.path.exists(out_path):
        drv.Delete(out_path)
    out_rast = drv.Create(out_path,xsize,ysize,1,gdal.GDT_Byte)
    out_rast.SetGeoTransform(geotrans)
    out_rast.SetProjection(rast_srs.ExportToWkt())

    # Rasterize filtered layer into the mask tif
    gdal.RasterizeLayer(out_rast, [1], mask_lyr, burn_values=[1],options=['ALL_TOUCHED=TRUE'])

    # Invert with Numpy - option 'INVERT=TRUE' had no effect
    out_rast.FlushCache()
    out_band = out_rast.GetRasterBand(1)
    out_arr = out_band.ReadAsArray()
    out_arr = np.logical_not(out_arr).astype(np.uint8)
    out_band.WriteArray(out_arr)
    out_band.ComputeStatistics(False)

    out_rast = None
    mask_ds = None
    rast_ds = None

The corresponding gdal_rasterize command is below:

gdal_rasterize GSHHS_f_L1.shp gdal_rasterize_mask.tif -burn 1 -at -i -l GSHHS_f_L1 -te 10.202 42.997 13.764 44.723 -ts 364 560`

Results

The output of gdal_rasterize is as expected with 0 values over land and 1 over the sea:

Output of gdal_rasterize

However, the output of my function only masks out the coastlines:

Output of gdal.RasterizeLayer

I've run the same function on a sample of other SAR scenes over different locations and it works fine, successfully masking all pixels that fall on land. The issue seems to be with this location.

I've ran this comparison using GDAL 2.1.2 and 2.2.3 with the same results.

Question

Is there a problem with my code, or has anyone noticed similar behaviour?

  • 1
    Strange error if it was working for other regions. If you want 1 for the ocean and 0 for land it would be interesting to see what happens if you fill your band with 1 out_rast.GetRasterBand(1).Fill(1) and then burn in the value 0 with rasterizing gdal.RasterizeLayer(out_rast, [1], mask_lyr, burn_values=[0],options=['ALL_TOUCHED=TRUE']) – Jascha Muller Mar 13 at 14:11
  • @JaschaMuller I tested again with this approach and got the same output from my function. Nice suggestion with the band.Fill(1) beforehand, i'd forgotten that was an available band method. – Ali Mar 13 at 15:19
  • Bummer, qualitatively it looks like the shape layer is picked up as a line rather than a polygon then. If you can give me the vital stats of your raster (geotransform, xsize, ysize and epsg) I can have a look in more detail. Maybe I can pick something up. – Jascha Muller Mar 13 at 17:20
  • Right! I had previously checked the geometries of the features after applying the spatial filter and they are all polygons. I've added the source raster details now so would be interesting to know if you can replicate this behaviour. – Ali Mar 14 at 20:57
0

I tried an alternative approach where, rather than setting a spatial filter on my mask layer and rasterising that result, I decided to use ds.ExecuteSQL on my mask data source object.

I constructed a spatialite query to find the intersection of my mask_ds geometries with the bounding box of my rast_ds. I then used this result layer to pass to the gdal.Rasterize function which successfully produced the land mask.

This doesnt quite answer the "why" but it looks like the problem was with the lyr.SetSpatialfilter method.

I've provided my working script below and included improvements such as the use of gdal.ApplyGeoTransform to get the corner coordinates of my raster as well as filling my output raster band with 1 and burning 0 (as suggested by Jascha in the comments).

import os
import numpy as np
from osgeo import gdal, ogr, osr
from timeit import default_timer as timer

def create_vector_mask(src_rast,mask_vect,dst_path,vect_lyr=0):
    '''
    Generates a mask raster for src_rast using the vector geometries of mask_vect.

    The mask is generated to mask out the interior of the vector geometries
    intersecting the raster i.e. 0 for areas within the geometry e.g. land
    cover and 1 for areas outside of the geometry e.g. sea.

    TODO: Add option to mask the exterior of the geometries.
    TODO: Add check for matching EPSGs and simplify query.

    Args:
        src_rast(str): Path to source raster file.
        mask_vect(str): Path to mask vector file.
        dst_path(str): Path to destination file.
        vect_lyr(int or str): Layer index or name to use for masking.
            Default is 0, the first layer.

    Returns:
        None
    '''

    # Get raster metadata
    rast_ds = gdal.Open(src_rast)
    if rast_ds is None:
        sys.exit('Failed to open src_file "{}"'.format(src_rast))
    gt = rast_ds.GetGeoTransform()
    b = rast_ds.GetRasterBand(1)

    # Get vector metadata
    mask_ds = ogr.Open(mask_vect)
    if mask_ds is None:
        sys.exit('Failed to open mask_vect "{}"'.format(mask_vect))
    mask_lyr = mask_ds.GetLayer(vect_lyr)
    if mask_lyr is None:
        sys.exit('Failed to open mask_vect layer "{}"'.format(vect_lyr))

    # Get EPSG info
    rast_srs = osr.SpatialReference(wkt=rast_ds.GetProjection())
    rast_srs.AutoIdentifyEPSG()
    rast_epsg = rast_srs.GetAttrValue('AUTHORITY',1)

    mask_srs = mask_lyr.GetSpatialRef()
    mask_srs.AutoIdentifyEPSG()
    mask_epsg = mask_srs.GetAttrValue('AUTHORITY',1)

    # Get raster corner points
    ul = gdal.ApplyGeoTransform(gt,0,0)
    ur = gdal.ApplyGeoTransform(gt,b.XSize,0)
    lr = gdal.ApplyGeoTransform(gt,b.XSize,b.YSize)
    ll = gdal.ApplyGeoTransform(gt,0,b.YSize)

    # Generate and execute spatial query on mask_ds
    rast_envelope = '\'POLYGON(({}))\''\
        .format(', '.join(['{} {}'.format(x,y) for x,y in [ul,ur,lr,ll,ul]]))

    sql_select = 'select st_transform(st_intersection(geometry, st_transform' \
        '(st_geomfromtext({0},{1}),{2})),{1}) as geometry'\
        .format(rast_envelope,rast_epsg,mask_epsg)

    sql_where = 'st_intersects(geometry,st_transform(st_geomfromtext({0},{1})'\
        ',{2}))'.format(rast_envelope,rast_epsg,mask_epsg)

    res_lyr = mask_ds.ExecuteSQL('{} from GSHHS_f_L1 where {}'\
        .format(sql_select,sql_where),dialect='sqlite')

    # Create raster to store mask
    drv = gdal.GetDriverByName('GTiff')
    if os.path.exists(dst_path):
        drv.Delete(dst_path)
    mask_rast = drv.Create(dst_path,b.XSize,b.YSize,1,gdal.GDT_Byte,
                        options=['TILED=YES','COMPRESS=DEFLATE','NBITS=1'])
    mask_rast.SetGeoTransform(gt)
    mask_rast.SetProjection(rast_srs.ExportToWkt())
    mask_band = mask_rast.GetRasterBand(1)
    mask_band.Fill(1)

    # Rasterize filtered layer into the mask tif
    gdal.RasterizeLayer(mask_rast, [1], res_lyr, burn_values=[0],
                        options=['ALL_TOUCHED=TRUE'])
    mask_rast.FlushCache()
    mask_rast = None

    mask_ds.ReleaseResultSet(res_lyr)
    mask_ds = None
    rast_ds = None

if __name__ == '__main__':
    src = r"source_raster.tif"
    msk = r"GSHHS_f_L1.shp"
    dst = '_mask'.join(os.path.splitext(src))
    t1 = timer()
    create_vector_mask(src,msk,dst)
    t2 = timer()
    print('Mask generated in {:0.3f} seconds'.format(t2-t1))

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.