0

I am trying to use polygon coordinates to plot polygons on GeoTIFFs of satellite images using gdal or rasterio, yet somehow I keep failing. There is two things I want to do:

  1. create a mask with only the polygon as data and the rest as nodata, and
  2. overlay the mask on the image to see how accurate the polygons are.

I simply do not know how to create a mask and overlay it in gdal in python when starting from WKT or GeoJSON strings.

I have the GeoTIFF loaded (and can plot it), and have the polygon coordinates in wkt and geojson format (and also as a geometry object). I can also transform the coordinates to and from desired projections. But I simply have no idea how I can overlay them on an image in GDAL or Rasterio to check the accuracy of the polygon (or better yet: create a separate mask which I can plot on top of the image).

Here is an example of a polygon I have in wkt in EPSG 4326:

'POLYGON ((5.5937209 52.24012314 0,5.5936411 52.24048512 0,5.59413417 52.2415213 0,5.59434149 52.24151408 0,5.59463832 52.24015416 0,5.5937209 52.24012314 0))'

I've seen some potential answers using .shp files, but haven't managed to recreate them from my starting point. It is not impossible that the answer is already somewhere on the forums, but I've been reading and searching for a long time and can't find the right one.

I do know how to do it in QGIS using the quickWKT plugin which is very straightforward.

enter image description here

Eventually however, I want to create a process in python to use these polygons to create masks for smaller tiles of GeoTIFFs for modelling. I wouldn't mind going straight to creating masks using polygons and simply plot them with an alpha over the image.

Improve this question
3
  • Your question is a bit vague on what you're trying to end up with. Perhaps you're running into a coordinate system problem (eg your raster in Google Mercator and your polygon in WGS84/Geographic) or is it that you want to create a map or serve a page of the polygon over the raster and can't work out how to overlay to create a result?
    – Michael Stimson
    Commented Jan 8, 2020 at 3:31
  • Seems to me that you want to rasterize the polygon. You can burn a raster map with a fixed value within the polygon, and a nodata value outside it. Then, you can combine both rasters (the original true color image with the rasterized polygon) into just one raster. I'm pretty sure that you can do it with GDAL Python binding. But you must define excactly what you want to do.
    – Gabriel De Luca
    Commented Jan 8, 2020 at 5:17
  • Thanks both for the reply. There is two things I want to do: 1)create a mask with only the polygon as data and the rest as nodata 2) but I also overlay the mask on the image to see how accurate the polygons are. It is not a coordinate system problem: I simply do not know how to create a mask and overlay it in gdal in python when starting from WKT or GeoJSON strings.
    – TheBackwardsman
    Commented Jan 8, 2020 at 8:32

1 Answer 1

Reset to default
1

If you are interested in masking the data, I would recommend rioxarray.

An example of doing so can be found here.

Here is a targeted example for what you probably want to do:

from shapely.geometry import mapping
from shapely.wkt import loads
import rioxarray

geom = mapping(loads('POLYGON ((5.5937209 52.24012314 0,5.5936411 52.24048512 0,5.59413417 52.2415213 0,5.59434149 52.24151408 0,5.59463832 52.24015416 0,5.5937209 52.24012314 0))'))

rds = rioxarray.open_rasterio("path_to_file.tif", parse_coordinates=False)
masked = rds.rio.clip([geom], "EPSG:4326", drop=False)
masked.rio.to_raster("masked_file.tif")
Improve this answer
4
  • Thanks for the reply, this is useful! With some adaptations I've been able to use your code to create a 'masked' object. The masked object this code returns however is considerably smaller (rds = (3, 34090, 32848), masked = (3, 311, 136)). This is because what we've just done is cut the polygon out of the image (when i plot the masked object it shows only a very small image that is the polygon contents). Do you know of a way that creates a mask the size of the original GeoTIFF with nodata everywhere accept for the polygons I would feed this code? Im looking for it but unsure where to look
    – TheBackwardsman
    Commented Jan 8, 2020 at 21:00
  • Yes - just add drop=False to the clip command.
    – snowman2
    Commented Jan 8, 2020 at 21:43
  • Thanks again! This might work. Sadly, I am getting a memory error (my GeoTIFFs are quite large, though I didn't expect a mem error since a mask can't be that large), and have been trying to get around that but haven't succeeded yet. I'll keep trying and accept your answer (and post my code snippet) when I manage to get it to succeed and see it work. Since I want to do this for many images; if you know of an efficient way out of the top of your head to do this without a memory error I'd love to hear.
    – TheBackwardsman
    Commented Jan 8, 2020 at 22:42
  • Just edited the answer to use parse_coordinates=False when opening the file. That might help a bit with memory.
    – snowman2
    Commented Jan 9, 2020 at 0:43

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