Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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!

share|improve this question
1  
I think you can make one mask, and use it for all of the rasters. – kyle 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 8:13

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)

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.