4

Suppose I have a satellite image formatted as a .tif file and a .kml file that contains a series of polygons that define an area of interest in this image. Now I can run:

gdalwarp satellite.tif satellite_aoi.tif -cutline aoi.kml

This will extract a new .tif image that only has the area of interest specified by the .kml file.

Now instead, what I want to do is perform the same area cut out but for the output, I want the areas of interest to be a mask. That is, I want the output of those regions to be a constant white color and preferably an image with a single channel. Is it possible to do this with gdalwarp? I looked at the documentation and I could not find any command to do what I need.

Improve this question
3
  • It sounds like you want to rasterize your .kml file. I believe you can do that with gdal_rasterize, but if you want the extents to match the satellite.tif file, you'll have to set them. Otherwise the output will have the same extents as the set of polygons in the kml file.
    – Jon
    Commented Nov 2, 2017 at 16:18
  • Yeah I need the extents to match the satellite.tif file. I tried gdal_rasterize but wasn't sure how to get the extents to match.
    – jlcv
    Commented Nov 2, 2017 at 21:39
  • It looks like you've got a solution that works, so no need to bang your head on it. I would get the satellite raster extents in python using the geotransform and raster size. You can also set your resolution to match using the geotransform output. Finally, I would use the -tap switch when calling gdal_rasterize to make sure the pixels aligned. All in all, it would be considerably more code than what you've posted.
    – Jon
    Commented Nov 2, 2017 at 22:30

2 Answers 2

Reset to default
3

As far I know, it's not possible to accomplish the task you need only with gdalwarp. Here's a working solution in Python (it requires GDAL 2.1.0+):

from osgeo import gdal

vector_layer = "mask.shp"
raster_layer = "raster.tif"
target_layer = "mask.tif"

# open the raster layer and get its relevant properties
raster_ds = gdal.Open(raster_layer, gdal.GA_ReadOnly)
xSize = raster_ds.RasterXSize
ySize = raster_ds.RasterYSize
geotransform = raster_ds.GetGeoTransform()
projection = raster_ds.GetProjection()

# create the target layer (1 band)
driver = gdal.GetDriverByName('GTiff')
target_ds = driver.Create(target_layer, xSize, ySize, bands = 1, eType = gdal.GDT_Byte, options = ["COMPRESS=DEFLATE"])
target_ds.SetGeoTransform(geotransform)
target_ds.SetProjection(projection)

# rasterize the vector layer into the target one
ds = gdal.Rasterize(target_ds, vector_layer, burnValues = [1])

target_ds = None
Improve this answer
4
  • Is it possible to rasterize each individual polygon as a different burn value?
    – jlcv
    Commented Nov 4, 2017 at 22:00
  • So I am using an older version of gdal 1100100 and there doesn't seem to be gdal.Rasterize. There is gdal.RasterizeLayer but I cannot get it working. Does it still work for this version? I tried gdal.RasterizeLayer(target_ds, [1], target_layer, burn_values=[1]) and I get TypeError: in method 'RasterizeLayer', argument 4 of type 'OGRLayerShadow *'
    – jlcv
    Commented Nov 4, 2017 at 22:41
  • Alright I tried: ds = gdal.RasterizeLayer(target_ds, [1], ogr.Open(vector_layer).GetLayer(0), burn_values=[1]); target_ds = None and I get Process finished with exit code 139 (interrupted by signal 11: SIGSEGV)
    – jlcv
    Commented Nov 4, 2017 at 23:30
  • But yeah I have to use version 1.10.1 and can't really change it :S
    – jlcv
    Commented Nov 5, 2017 at 0:12
1

This is the best solution I can think of in Python (not the most efficient but works)...

shutil.copy('satellite.tif', 'ones.tif')
ds = gdal.Open('ones.tif', gdal.GA_Update)
data = ds.ReadAsArray()
im_mask = np.ones((data.shape[1:])) * 255
ds.GetRasterBand(1).WriteArray(im_mask)
ds.GetRasterBand(2).WriteArray(im_mask)
ds.GetRasterBand(3).WriteArray(im_mask)
ds = None

command = 'gdalwarp ones.tif satellite_aoi.tif -cutline aoi.kml'
os.system(command)
Improve this answer
1
  • Although for this solution I output 3 bands when I really only need one...
    – jlcv
    Commented Nov 2, 2017 at 16:02

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.