1

I have a raster, and I'm trying to create a polygon and export to a shapefile that is showing one solid layer between lets say 1000m and 3000m in elevation.

So it seems like I need to filter out all the data outside of those elevations, but I'm not sure how to go about that.

How do I do this using GDAL in Python?

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
m0ngr31
  • 211
  • 1
  • 7
  • 3
    You could try gdal_calc -A your_raster -calc"logical_or(A>1000,A<3000) " based on this post https://gis.stackexchange.com/questions/204387/gdal-calc-raster-calculator-syntax-for-logical-operators-and-other-functions that should make a binary raster with 0 and 1, polygonize the binary raster - it's less intensive than a continuous raster. – Michael Stimson May 09 '18 at 04:11
  • What does the A represent? – m0ngr31 May 09 '18 at 04:52
  • 1
    A is your raster, specified by -A.. this allows you to have multiple rasters in the calculation -A first_raster -B second_raster etc.. then in the calc you can specify A operator B (-calc "A+B" for example), read about it http://www.gdal.org/gdal_calc.html – Michael Stimson May 09 '18 at 05:19
  • Okay that makes sense. It appears I'm doing something wrong with gdal_calc, because it fills the whole state as black. If I do it in QGIS with Raster Calculator idaho@1 >= 1000 AND idaho@1 <3000 that seems to create the right thing. But when I got to rasterize it to shapefile, it just does a line instead of a polygon. I'm sure I'm missing something stupid. – m0ngr31 May 09 '18 at 05:46
  • 1
    Is your gdal raster binary? How are you exporting to a shapefile? This might help https://docs.qgis.org/2.8/en/docs/training_manual/complete_analysis/raster_to_vector.html get you back on track. – Michael Stimson May 09 '18 at 05:59
  • I don't know what happened, but it's working now. Thank you so much! You're a live saver. If you want to add an actual answer I can make that the accepted answer. – m0ngr31 May 09 '18 at 06:07
  • 2
    Feel free to put in your own answer with screen shots and links based on any direction you've gleaned from my comments and I'll upvote it tomorrow when I get back in. – Michael Stimson May 09 '18 at 06:09

1 Answers1

5

After playing around with the links Michael gave me above, this ended up being the solution:

gdal_calc.py -A idaho.tif --outfile=idaho-filtered.tif --calc="1*logical_and(A>=1000,A<3000)" --NoDataValue=0
gdal_polygonize.py idaho-filtered.tif -f "ESRI Shapefile" idaho-filtered.shp idaho-filtered

Here's the results...

Original: enter image description here

After gdal_calc: enter image description here

Shapefile: enter image description here

m0ngr31
  • 211
  • 1
  • 7
  • 1
    Excellent, especially the section about --NoDataValue=0 to eliminate the polygons you're not interested in. If you have the time can you put in a screen shot of the original idaho.tif, binary idaho-filtered.tif and the resulting polygon. – Michael Stimson May 09 '18 at 22:23
  • Thank you! I had a lot of trouble finding how to do that! – Tristram Gräbener Dec 30 '19 at 18:18