A .tif image has been created using the NDVI formula in QGIS. The resulting image shows a region, with varying grades of NDVI values over the image. A .shp file has then been created for an area of interest, then the corresponding _mask file, using rasterio. I've been unable to identify the necessary steps in Python to firstly, remove all values on the edge (i.e. overlapping) the polygon's edge and to remove all "dead" values (dummy, None, whatever), then secondly, I need to look at every pixel in this _mask file and identify the NDVI value associated with it. Then sum all values and average it. My code so far:

def __init__(self):
    self.shape_file = r"path_1"
    self.raster_file = r"path_2"
    self.mask_file = r"path_3"

    out_image, no_data = self.create_mask_image()
    self.summarise_raster_values(out_image, no_data)

def create_mask_image(self):
    with fiona.open(self.shape_file, "r") as shapefile:
        features = [feature["geometry"] for feature in shapefile]

    # Extract raster values values within the polygon
    with rasterio.open(self.raster_file) as src:
        if src is None:
            print ("Raster file could not be opened")
        out_image, out_transform = rasterio.mask.mask(src, features, crop=True)
        out_meta = src.meta.copy()

        # BS test code
        x = (src.bounds.left + src.bounds.right) / 2.0
        y = (src.bounds.bottom + src.bounds.top) / 2.0
        for val in src.sample([(x, y)]):
            print (val)

    # No data values of the original raster (optional)
    no_data = src.nodata
    print (no_data)

    out_meta.update({"driver": "GTiff",
                     "height": out_image.shape[1],
                     "width": out_image.shape[2],
                     "transform": out_transform})
    with rasterio.open(self.mask_file, "w", **out_meta) as dest:

    return out_image, no_data

def summarise_raster_values(self, out_image, no_data):
    # Extract the values of the masked array
    data = out_image.data[0]
    # Extract the row, columns of the valid values
    row, col = np.where(data != no_data)
    #elev = np.extract(data != no_data, data)

The problem is, how to make the calculations within my summarise_raster_values method. Any ideas how I can achieve this?

The output of print (data):

[[-3.4028235e+38 -3.4028235e+38 -3.4028235e+38 ... -3.4028235e+38
  -3.4028235e+38 -3.4028235e+38]
 [-3.4028235e+38 -3.4028235e+38 -3.4028235e+38 ... -3.4028235e+38
  -3.4028235e+38 -3.4028235e+38]
 [-3.4028235e+38 -3.4028235e+38 -3.4028235e+38 ... -3.4028235e+38
  -3.4028235e+38 -3.4028235e+38]

I was expecting NDVI values here, but as you can see the values are exceedingly large.

  • why not look at existing zonal statistic tools? docs.qgis.org/2.18/en/docs/user_manual/plugins/… – radouxju Nov 9 '18 at 15:23
  • If you want to use Python, just install rasterstats. It wraps around rasterio and does exactly what you want, including returning masked arrays. – Jon Nov 9 '18 at 15:39
  • @Jon: that's a great suggestion, I think that is exactly what I was looking for! I just checked the code and the results from the stats. It appears the min value is a negative number. When I upload the same NDVI_mask file into QGIS, the values are positive (about 0.1 to 0.7). However rasterstats produces: [{'max': 0.7470259070396423, 'sum': 1222844.25, 'median': 0.46611133217811584, 'majority': 0.5, 'min': -0.14378099143505096}]. I'm wondering how this negative value is occurring. – pymat Nov 9 '18 at 16:13
  • Is your min value the nodata value? If so, just pass it in to rasterstats explicitly so it ignores it. I've done checks against QGIS Zonal Stats and rasterstats and they're in perfect agreement, except when there are only a few pixels in the polygon of interest. – Jon Nov 9 '18 at 19:17
  • @Jon: my src.nodata comes out at -3.40282346639e+38, which is different to the nodata value rasterstats produces. – pymat Nov 10 '18 at 9:06

Your Answer

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

Browse other questions tagged or ask your own question.