1

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.

  • 1
    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
2

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[0],
        'width': data.shape[1],
        '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

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.