# Comparing pixel values at a given point in two rasters

I have two raster images, one for land use and one for NDVI, i would like to calculate LAI based on land use. So i need to compare the pixel value in the NDVI raster and Land use raster then select which formula to use. How do i approach this in python? I am a beginner.

• What do you want as an output? Do you just need to compare the two values somewhere, or do you want to calculate a new raster? – Marc Pfister Dec 11 '18 at 21:57
• Do you have any code to start from? Do both rasters have the same cell size, extent and origin? Can you give an example of a couple of NDVI and Land Use combinations and their related formula? It would make it easier to give example code if there were a few conditions taken to conclusion. – Michael Stimson Dec 11 '18 at 22:45

## 1 Answer

The question you are asking is quite complex to answer comprehensively, but I can give you some structure and you can fill the gaps. My answer is using `arpcy`, that comes with ArcGIS. You can achieve something very similar with `rasterio` (see commented lines)

First of all, I always like to work directly with the numpy arrays. This can be achieved in arcpy by:

``````import arcpy

raster_path = 'C:/blablabla/raster_bla.tif'
raster = arcpy.Raster(raster_path)
raster_data = arcpy.RasterToNumPyArray(raster)
# import rasterio
# raster = rasterio.open(raster_path)
# raster_data = raster.read()
``````

raster_data is now a numpy.ndarray containing all the pixel values. Before proceeding to the next steps, you should perform some tasks to make sure everything is what you expect:

• load both rasters (of course!)
• using the `meta = arcpy.Describe(raster)` method you can check things like width, height, pixel_size and others to make sure that both rasters are actually aligned perfectly. If this is not the case things will not work. Some background on the Describe method can be found here
• You may want to check that your arrays are 2D (e.g. the have only one band)

Next step would be something like:

``````def leaf_area_index(land_use, ndvi):
result = ndvi.copy()
for lu in np.unique(land_use):
results[land_use == lu] = _calculate_lai(ndvi[land_use == lu], lu)
return result
``````

where:

• `land_use == lu` is pretty much selecting pixels belonging to each of your land use values
• `_calculate_lai()` should be a function that gets the NDVI values and the land use code, and returns the leaf_area_index (I suspect each land use code will have different exponents to relate ndvi to lai)

After that you probably want to save the output raster. In arcpy this can be naively achieved by using something like this:

``````def np_to_raster(array, raster_ref):
"""Utility to convert a numpy array to a Raster object using another Raster as reference

Arguments:
-----------
array : np.ndarray
data to convert
raster_ref : arcpy.Raster, or str
reference raster or path to reference raster

Returns:
-----------
out : arcpy.Raster
converted Raster object
"""
if not isinstance(array, np.ndarray):
raise TypeError('The input must be a numpy array')
meta = arcpy.Describe(raster_ref)
cell_size_x = meta.meanCellWidth
cell_size_y = meta.meanCellHeight
lowerleftcorner = meta.extent.lowerLeft

out_raster = arcpy.NumPyArrayToRaster(array, lowerleftcorner, cell_size_x, cell_size_y, 0)
return out_raster
``````

this should be called by your main function to covert the numpy.array to arcpy.Raster object.

in rasterio you can do

``````def write_img(array, raster_ref, outpath):
newprofile = {
'count': 1,
'crs': raster_ref.crs,
'dtype': 'uint8',
'driver': 'GTiff',
'transform': raster_ref.transform,
'height': data.shape,
'width': data.shape,
'blockxsize': 256,
'tiled': True,
'blockysize': 256,
'nodata': None
}

with rasterio.open(outpath, 'w', **newprofile) as dst:
dst.write(data)
``````

In summary (this is NOT a working code, you need to fill the gaps)

``````def run(land_use_path, ndvi_path, output_path):

#here you load your stuff
land_use = arcpy.Raster(land_use_path)
lu_meta = arcpy.Describe(land_use)

#here you test the rasters to see if they are aligned
assert lu_meta.height == ndvi_meta.height

#here you get the data out of the raster
lu_data = arcpy.RasterToNumPyArray(land_use)

#here you do your calculation
lai_data = leaf_area_index(lu_data, ndvi_data)

#here you convert the result to raster
lai_raster = np_to_raster(lai_data, land_use)
lai_raster.save(output_path)
# or write_img(lai_data, land_use, output_path) if you are working with rasterio
``````

Hope this will help!

• That's all very nice but the OP doesn't appear to have Esri software, the question is about GDAL in python. – Michael Stimson Dec 12 '18 at 0:20
• rasterio is pretty much gdal – Alessio Arena Dec 12 '18 at 0:29