1

Does anyone have a method for extracting pixels from a GeoTIFF that are of a particular value. So for example if a band had values from 0 to 10 and I wanted to make an image just using the zero values and to leave the other pixel values as unknown.

My GeoTIFF is 19 GB large so I'd prefer to avoid using QGIS. I've already tried GRASS and it would freeze on me. I would prefer a solution that can use the GDAL tools but Python would be good too.

My eventual goal is to extract the polygons that surround regions of interest using the monochrome image generated from the above.

  • Have you tried raster calculator? – BERA May 6 at 12:19
  • @BERA I have but it was freezing. – Janusz May 6 at 12:20
  • Or have you tried gdal_calc python script gdal.org/programs/gdal_calc.html? – user30184 May 6 at 12:20
  • @user30184 no I haven't. I will try that right now. – Janusz May 6 at 12:22
  • 1
    It will be faster to test with an image that is smaller than 19 GB. You can find the limits later and with so big image you will want to use a few creation options as well, at least tiled=yes and some compression method. – user30184 May 6 at 12:24
3

I've not tried this but I think the example in the gdal_calc manual page

gdal_calc.py -A input.tif --outfile=result.tif --calc="A*(A>0)" --NoDataValue=0

could be changed to

gdal_calc.py -A input.tif --outfile=result.tif --calc="A*(A==1)" --NoDataValue=0

should work to extract value 1.

| improve this answer | |
  • I'm running it, we'll see if I get the desired result. – Janusz May 6 at 12:26
  • That worked for what I needed. In particular I used --calc="A==0" since I only wanted to extract the border around the data. – Janusz May 6 at 14:20
1

this is one way to do what you want using GDAL and NUMPY.

from osgeo import gdal
import numpy as np

rasterFile = 'C:/Users/pepito/yourinputraster.tif'
in_ds = gdal.Open(rasterFile)
in_band = in_ds.GetRasterBand(1)
in_data = in_band.ReadAsArray()

#this is the line that makes the reclass
reclas = np.where(in_data == 0, 1, 0)

# output as TIFF
gtiff_driver = gdal.GetDriverByName('GTiff')
out_ds = gtiff_driver.Create('C:/Users/pepito/youroutputraster.tif',in_band.XSize, in_band.YSize, 1, in_band.DataType)
out_ds.SetProjection(in_ds.GetProjection())
out_ds.SetGeoTransform(in_ds.GetGeoTransform())
out_band = out_ds.GetRasterBand(1)
out_band.SetNoDataValue(0)
out_band.WriteArray(reclas)
out_ds.FlushCache()

However if your raster is 19GB you will probably have to process it in chunks.

| improve this answer | |

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.