4

How do I exclude no data values when doing a raster calculation with rasterio. I realise that the files are read into an array, and that "0" is a nodata value. However, when I do a simple sum like "image + 1" it applies the math to the nodata value as well. So what is the best way to 'ignore' nodata values in the calculation?

FYI: I am calibrating the Landsat bands to reflectance and this is what my current snippet of the band math looks like:

    if "B1.TIF" in band:
        print band
        image = "%s" % (band)
        with rasterio.open(image) as src:
            image_read = src.read(1)
        #np.seterr(divide='ignore', invalid='ignore')
        constant = 0.01745329251994444444444444444444 #Constant is calculated (3.14/180) which is converting the sun-angle to sun_radians which was suggested by WOlfgang
        set_mb = "{:.16f}".format(float(mult_band1))
        set_ab = "{:.16f}".format(float(add_band1))
        toa_1 = (float(set_mb) * image_read.astype(float)) + (float(set_ab))
        solar_z = np.cos((90-float(sun_elevation))*float(constant))
        toa_2 = (toa_1.astype(float) / solar_z) * 100000
        kwargs = src.meta
        with rasterio.open(outfile, 'w', **kwargs) as dst:
            dst.write_band(1, toa_2.astype(rasterio.uint16))
        #print image_read.mask, toa_2

Is there a simple way to ignore nodata values in the calculation?

3

You can get a Numpy masked array that covers up nodata values from Rasterio by adding a keyword argument: src.read(1, masked=True). Operations on a masked array do not use the covered up elements.

If your dataset has no defined nodata value, but you want to use for example 0, read out a non-masked array and mask it yourself:

image_read = src.read(1)
image_read_masked = numpy.ma.masked_array(image_read, mask=(image_read == 0))
  • thanks for your response. It seems like that the band math is still done on the 0 values when adding that keyword argument. 'src = rasterio.open(image) image_read = src.read(1, masked=True) image_read.mask print image_read' '[[0 0 0 ..., 0 0 0] [0 0 0 ..., 0 0 0] [0 0 0 ..., 0 0 0] ..., [0 0 0 ..., 0 0 0] [0 0 0 ..., 0 0 0] [0 0 0 ..., 0 0 0]]' – Jens Hiestermann Jan 11 '17 at 7:59
  • ok so I had a look at the image nodatavals using src.nodatavals and it returned (None,) And that is probably why src.read(1, masked=None) is not working as a mask. Any idea as to how to update the bands nodataval to zero? – Jens Hiestermann Jan 11 '17 at 8:35
  • Right: if the GeoTIFF has no nodata value set, masked=True has no effect. In your case, instead of modifying the Landsat file's nodata value do it all in code: image_read = src.read(1) and then image_read_masked = numpy.ma.masked_array(image_read, mask=(image_read == 0)). – sgillies Jan 11 '17 at 10:54
  • 1
    Another way is with numpy.ma.masked_equal(x, value) – Mike T Jan 11 '17 at 22:22
  • @sgillies Could I change the nodata value from "None" to "-9999"? ` then image_read = np.ma.masked_values(band, -9999)` would have the desired result. I need an eloquent way to mask out nodata value which is 'None' from my raster calculations. Please help – Jens Hiestermann Jan 12 '17 at 6:53

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.