Terrain correction

Functions for implementing terrain correction algorithm (credit to Wim Nursal and team, LAPAN)

Original code at https://github.com/Forests2020-Indonesia/Topographic-Correction/blob/master/Topographic%20Correction.py Method b: https://ieeexplore.ieee.org/document/8356797

pyeo.terrain_correction.build_sample_array(raster_array, slope_array, red_band_index, ir_band_index)

Returns a set of pixels in raster with slope > 18deg + ndvi > 0.5

pyeo.terrain_correction.calculate_ic_array(slope_raster_path, aspect_raster_path, raster_datetime, ic_raster_out_path=None)

Given a slope and an aspect raster, creates an array of the illumination conditions as specified in https://ieeexplore.ieee.org/document/8356797, equation 9. The Pysolar library is used to calculate solar position.

Parameters
  • dem_raster_path – The path to a raster containing the DEM in question

  • raster_datetime – The time of day _with timezone set_ for the

  • ic_raster_out_path – If present, saves a raster of the illumination condition

Returns

  • An array of illuminatiton conditions for every pixel in the input DEM.

  • Each pixel is a value between -1 and 1

pyeo.terrain_correction.do_terrain_correction(raster_path, dem_path, out_raster_path, raster_datetime, is_landsat=False)

Corrects for shadow effects due to terrain features. Algorithm:

  • Generate slope and aspect from DEM using gdaldem

  • Calculate solar position from datatake sensing start and location of image

  • Calculate the correction factor for that image from the sun zenith angle, azimuth angle, DEM aspect and DEM slope

  • Build a mask of green areas using NDVI

  • Perform a linear regression based on that IC calculation and the contents of the L2 image to get ground slope(?)

  • Correct pixel p in original image with following: p_out = p_in - (ground_slope*(IC-cos(sun_zenith)))

  • Write to output

Parameters
  • raster_path – A path to the raster to correct.

  • dem_path – The path to the DEM

  • out_raster_path – The path to the output.

  • raster_datetime – A datetime.DateTime object with timezone set

pyeo.terrain_correction.generate_slope_and_aspect_rasters(dem_raster_path, out_directory)

Using GDAL, splits a DEM into slope and aspect rasters

pyeo.terrain_correction.get_dem_slope_and_angle(dem_path, slope_out_path, aspect_out_path)

Produces two .tifs, slope and aspect, from an imput DEM. Assumes that the DEM is in meters.

pyeo.terrain_correction.get_pixel_latlon(raster, x, y)

For a given pixel in raster, gets the lat-lon value in EPSG 4326.