2

I think this should be a relatively simple process, but I'm not quite sure how to approach it.

I have a shapefile containing 1000's of polygons. I also have 19 rasters. All of the polygons fall withing the extent of the rasters. The rasters are equal in extent and resolution. Previously I have used GDAL to extract data for each of these polygons from each of the 19 rasters. This worked very well. However, I now wish to extract the pixels into xyz format or similar for all the areas not covered by these polygons.

Currently I am using gdal_rasterize to burn a chosen nodata value into each of the rasters indicating where the polygons occur. I would then need to translate the raster removing all of the nodata values. This is a very heavy process and is taking an awefully long time for 1 raster let alone 19.

My question is - is there a way to gdal_rasterize all of my polygons onto just 1 of the rasters, then use the nodata values from this raster as a mask, and apply these nodata values to the same locations on the other 18 rasters. I'm sure this would be quicker than my current method, but I don't know how to do this in GDAL.

Any help will be very much appreciated!

JPD
  • 1,678
  • 2
  • 19
  • 33
  • 1
    I think you can make one mask, and use it for all of the rasters. –  Jan 30 '13 at 16:07
  • What would the command be to apply the nodata values from 1 raster to the other then? – JPD Jan 31 '13 at 08:13

1 Answers1

1

I have a shapefile containing 1000's of polygons. I also have 19 rasters. All of the polygons fall within the extent of the rasters. The rasters are equal in extent and resolution. Previously I have used GDAL to extract data for each of these polygons from each of the 19 rasters. This worked very well. However, I now wish to extract the pixels into xyz format or similar for all the areas not covered by these polygons

Make a single polygon (sp) that covers the extent of all polygons (p1, p2) and the rasters (which are all the same extent so it will work). Then create a model or similar.

  1. Erase sp by p1 = p11 iterate for all polygons so p2 = p21, p3 - p31.
  2. Now you can just use the same method as before to extract your values for all the information outside (p11, p21, p31,etc)

This should be very rapid, indeed part 2 should take the same time as before assuming outside size = inside size and part one is a simple vector function so can be done rapidly.

Now it looks like you are on a open source path so just use Difference in QGIS as opposed to Erase in ArcGIS.

I am unclear in my mind how to do part 1 in OGR but I did find that the OGR function to Erase a geometry is named Difference. I found this Python API document,[osgeo.ogr Python] (http://gdal.org/python/osgeo.ogr-module.html)

Community
  • 1
If you do not know- just GIS
  • 7,416
  • 28
  • 71