Latest Contributions by PeterWilson

Hi Christine Dartiguenave‌ Thanks for the update, really appreciate you taking the time to get back to me. I'll wait a bit longer then thanks ... View more

Hi Rodrigo Rocha  I'll have a look at the following and get back to you shortly. Sorry I've not been on the forums either for a while. I've changed my career direction from geospatial to data science. Regards Peter ... View more

Hi ArcHydro Community  I've been using ArcGIS 10.3.1 for quite a while with ArcHydro 10.3.0.78 and due to my python package that I've developed have been reluctant to upgrade to ArcGIS 10.6. I wanted to find out from Christine Dartiguenave‌ if there are any major changes from ArcHydro 10.3 to 10.6 version that would break my existing code, before making the leap. ... View more

Thanks Dan , will check it out and get back to you. ... View more

Thanks Dan ... View more

Hi Dan I completed writing my Python module to automatically generate NDVI rasters from Landsat 8 tiles (Downloaded from USGS EarthExplorer) with TOA Reflectance and Sun Angle correction. NDVI: With TOA Reflectance and Sun Angle correction NDVI: Value range (between 1 and -1) ''' Created on 23 Sep 2017 Create NDVI Rasters with TOA Reflectance and Sun Angle correction @author: PeterW ''' # import site-packages and modules import re import argparse from pathlib import Path import arcpy from arcpy.sa import * def list_landsat_tiles(landsat_dir):     """     Create a list of Landsat 8     tiles bands 4 & 5.     """     # Determine how to prevent nested loops - Big 0 Notation     landsat_tiles = {}     ndvi_bands = ['_B4.TIF', '_B5.TIF', '_MTL.txt']     p = Path(landsat_dir)     for directory in p.iterdir():         # add error handling to validate its subdirectories         for band in ndvi_bands:             match = '*{0}'.format(band)             landsat_tile = str(directory.glob(match).next())             landsat_name = str(directory.stem)             landsat_key = re.findall('_(\d{6})_\d{8}_(\d{8})', landsat_name)[0]             landsat_tiles.setdefault(landsat_key, []).append(landsat_tile)     return landsat_tiles def remove_zero_values(landsat_tiles):     """     Convert zero cell values     to NoData.     """     arcpy.CheckOutExtension('Spatial')     for k, v in landsat_tiles.iteritems():         red_band = SetNull(v[0], v[0], 'value=0')         NIR_band = SetNull(v[1], v[1], 'value=0')         v[0] = red_band         v[1] = NIR_band     arcpy.CheckInExtension('Spatial') def extract_reflectance_coefficients(landsat_tiles):     """     Extract the reflectance     coefficients from metadata     txt file.     """     for k, v in landsat_tiles.iteritems():         with open(v[2]) as mlt:             lines = mlt.read().splitlines()             reflect_mult = float(lines[187].split('=')[1])             v[2] = reflect_mult             reflect_add = float(lines[196].split('=')[1])             v.append(reflect_add)             sun_elev = float(lines[76].split('=')[1])             v.append(sun_elev) def toa_reflectance_correction(landsat_tiles):     """     Correct landsat 8     bands 4 & 5     for TOA reflectance     """     arcpy.CheckOutExtension('Spatial')     for k, v in landsat_tiles.iteritems():         reflect4 = (v[2]*v[0])+v[3]         v[0] = reflect4         reflect5 = (v[2]*v[1])+v[3]         v[1] = reflect5     arcpy.CheckInExtension('Spatial') def sun_angle_correction(landsat_tiles):     """     Correct Landsat 8     bands 4 & 5     for sun angle     """     arcpy.CheckOutExtension('Spatial')     for k, v in landsat_tiles.iteritems():         sun4 = (v[0]/(Sin(v[4])))         v[0] = sun4         sun5 = (v[1]/(Sin(v[4])))         v[1] = sun5     arcpy.CheckInExtension('Spatial') def calculate_ndvi(landsat_tiles, output_dir):     """     Generate NDVI from     preprocessed     landsat 8 bands 4 & 5     """     arcpy.env.overwriteOutput = True     arcpy.CheckOutExtension('Spatial')     for f, v in landsat_tiles.iteritems():         NDVI_name = '_'.join(f)         arcpy.AddMessage('Processing {0}.tif NDVI'.format(NDVI_name))         Num = Float(v[1] - v[0])         Denom = Float(v[1] + v[0])         NDVI_raster = Divide(Num, Denom)         NDVI_output = '{0}\\{1}.tif'.format(output_dir, NDVI_name)         NDVI_raster.save(NDVI_output)     arcpy.CheckInExtension('Spatial') def main(landsat_dir, output_dir):     """     Determine NDVI for     each Landsat tile.     """     landsat_tiles = list_landsat_tiles(landsat_dir)     remove_zero_values(landsat_tiles)     extract_reflectance_coefficients(landsat_tiles)     toa_reflectance_correction(landsat_tiles)     sun_angle_correction(landsat_tiles)     calculate_ndvi(landsat_tiles, output_dir) if __name__ == '__main__':     parser = argparse.ArgumentParser(description='Calculate the NDVI for Landsat 8 tiles')     parser.add_argument('--landsat_dir', metavar='path', required=True,                         help='Input Landsat 8 tile directory')     parser.add_argument('--output_dir', metavar='path', required=True,                         help='Output NDVI directory')     args = parser.parse_args()     main(landsat_dir=args.landsat_dir,          output_dir=args.output_dir) ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ The only question that I have is that my output raster is between 1 and -1 but is not represented as 1 to -1 is that because of the actual values within my tile or is there a problem with my calculation? ... View more

Hi Cody Thanks for the response. I'd like to determine the NDVI for the dry and wet season for my study area and then determine the difference between the two. I'm currently still using ArcGIS 10.3.1, unfortunately no ArcGIS Pro. The link that I provided in my initial post is a blog post by ESRI Australia: Using Python to calculate NDVI with multiband imagery, that has a Python script for determing the NDVI. If you are able to assist in the approach I need to following once I've determined the NDVI for say (i.e. 2017 Feb (Wet) and 2017 Oct (Dry)) ,how I can show the difference in NDVI between the wet and dry season. I'd like to be able to code the following in Python as its part of a larger water project that I'm busy with. ... View more

Hi Dan I found the following on ResearchGate: Based on the following I presume that I should be able to calculate the NDVI using same formula, its just that the Landsat 8 is being stored in 16-bit and not 8-bit. ... View more

Hi Richard Thanks I clean forgot about Python's Pillow module. I was able to batch all the images and correct them. ... View more