0

I have a differential surface model dataset. That was made: DSM subtracted from DTM, so for example trees have negative height. It is a raster (TIFF) dataset, Type=Float32.

I want to create a shapefile that will have three different values:

  • 0: error (all the pixels in the original dataset with values > 0)
  • 1: OK (all the pixels in the original dataset with values < 0)
  • 2: nodata (all the pixels in the original dataset with values = -10000000000)

I assume first thing to create a raster with pixel values 0, 1, 2 only and than vectorize it. Probably the vectorization will not be a problem, using gdal_polygonize.py but creating the raster is quite a challenge for me, although I was able to achieved some minor success using gdal_calc.

I'd like to do this using GDAL, do you have an idea?

  • This type of operation is called a reclassification. There is a reclassification tool in ArcGIS but you would need a spatial analyst license. QGIS (free) has a Reclassify by table tool or you can use SAGA in QGIS to reclassify rasters. If you want to only use GDAL I found this link gis.stackexchange.com/questions/116473/…. Good luck. – GBG May 23 at 15:16
  • You can achieve this with gdal_calc and gdal_polygonize. If you want actual help, you should post what you've tried with gdal_calc and your results--i.e. describe what you mean by "minor success." If you use gdal within Python, you can just load the raster in as a numpy array and do the classification with numpy slicing, although you will need some extra coding to re-save as a geotiff. – Jon May 23 at 20:54
0

You can use the gdal python bindings to read the array and then numpy to change the pixel values. Then you can use gdal again to overwrite your original raster.

import gdal

fn = 'path_to_your_raster.tif'
ds = gdal.Open(fn, 1)   # open TIFF file in writing mode
arr = ds.ReadAsArray()  # read the data as a numpy array

# specify your conditions
con1 = (arr > 0)
con2 = (arr < 0)
con3 = (arr == -10000000000)

# change values according to the conditions
arr[con1] = 0
arr[con2] = 1
arr[con3] = 2

# overwrite the TIFF data
band = ds.GetRasterBand(1)
band.WriteArray(arr)

# save file and close
ds.FlushCache()
del ds, band
  • As far as I understand 'numpy' is installed correctly (based on this) but when I try to run your script I got several error messages: "ImportError: numpy.core.multiarray failed to import", "ImportError: No module named _gdal_array" – STO May 27 at 13:36
  • What Python version do you have? And is it a stand alone Python installation or are you using and installation that came with another programm (e.g. QGIS)? – Marcelo Villa May 27 at 14:06
  • v2.7 came with ArcGIS I think, but I have QGIS as well v 3.2.3 – STO May 28 at 9:36
  • @STO Can you try to run this script in the QGIS Python console? The Python installation that comes with ArcMap does not include GDAL. – Marcelo Villa May 28 at 16:40
1

If you can use R, it's raster package supports intuitive raster algebra that will make your task very simple. Try this:

library(raster)

rred <- brick(system.file("external/rlogo.grd", package="raster"))[['red']] ## MWE, use your data
values(rred) <- values(rred) - mean(values(rred))

rbad <- rred > 0
rgood <- rred <= 0
rmiss <- rred == -100 ## please use your true value of -10e9, this is only for MWE

writeRaster(rgood, filename='good.tif')
writeRaster(rbad, filename='bad.tif')
## etc

now polygonize them with GDAL, if you wish. Alternatively install and load gdalUtils and do it all inside R. The raster package can do polygonization too, but it's insufferably slow.

  • Sorry, I don't know what is "R". Can you give me a link? – STO May 24 at 9:31
  • Certainly. R was designed as a statistical computing language, and it's grown to be a capable GIS too. Download it here. It's package management is also superior to most--install the raster package with install.packages("raster") Load packages up like library(raster) – 0mn1 May 27 at 1:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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