POST
|
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
10-04-2018
09:03 AM
|
0
|
0
|
1635
|
POST
|
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
10-04-2018
07:07 AM
|
0
|
0
|
678
|
POST
|
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
10-04-2018
06:56 AM
|
0
|
1
|
1635
|
POST
|
I currently still using ArcGIS 10.3 and had to generate a 1m resolution DTM for the whole of Dar es Salaam for a Flood Risk Analysis task on one of our projects. The first problem is that the source of our input data was contour shapefiles at a 1m resolution. I initially created a 5m resolution DTM from the contour datasets for preliminary results which took almost 5 days to run. This was a problem as I had to get the resolution down to 1m to run my analysis. I generated a overlapping index grid allowing for 20% overlap and sliced the contour datasets into 500 tiles. I then generated a Python module utilising Python's Multiprocessing module to process the tiles across 16 cores on our servers. The DTM tiles were completed in under 12 hrs for the entire study area. The Python module which I've attached is below, is my first draft of utilizing Python's Multiprocessing. I'm currently re-writing it now that I have a better understanding how the Multiprocessing module works. NB: I'd like to thank Duncan Hornby for his post: Create a script tool that uses multiprocessing, that helped me create the following script Hope the following helps anyone else trying to achieve the same for generating large DTM's Dar es Salaam: 1m DTM Tiles '''
Created on 01 Feb 2018
Multi-threaded
Topo to Raster (3D)
@author: PeterW
'''
# import site-packages and modules
import multiprocessing
import argparse
import arcpy
from functools import partial
# set environment settings
arcpy.env.overwriteOutput = True
# check-out extensions
arcpy.CheckOutExtension('3D')
def dtm_work(tiles_directory, cont):
"""
Multiprocess
Topo to Rasster
for multiple
DTMs
"""
index_name = cont[0]
inContours = cont[1]
try:
output_dtm = '{0}\\{1}'.format(tiles_directory, index_name)
arcpy.TopoToRaster_3d(inContours,
out_surface_raster=output_dtm,
cell_size=1)
return True
except Exception as e:
print e.message
return False
def dtm_multi(index_grid, tiles_directory):
"""
Create overlapping
DTM tiles from
contour tiles &
index grid
"""
try:
cont_param = []
with arcpy.da.SearchCursor(index_grid, ['SHAPE@', 'PageName', 'Path']) as scur: # @UndefinedVariable
for row in scur:
index_name = row[1]
input_cont = '{0} {1} {2}'.format(row[2], 'ELEVATION', 'Contour')
cont_param.append([index_name, input_cont])
func = partial(dtm_work, tiles_directory)
arcpy.AddMessage('Sending to pool')
# declare the number of cores to use, use 2 less than the max
cpu_num = multiprocessing.cpu_count()-2
# create the pool object
pool = multiprocessing.Pool(processes=cpu_num)
# fire off list to worker function
# res is a list that is created with what ever the worker function returns
res = pool.map(func, cont_param)
pool.close()
pool.join()
# if an error has occurred report it
if False in res:
arcpy.AddError('A process\thread failed!')
arcpy.AddMessage('Finished Multiprocessing DTM Tiles')
except arcpy.ExecuteError:
# Geoprocessor through error
arcpy.AddError(arcpy.GetMessages(2))
if __name__ == '__main__':
parser = argparse.ArgumentParser(description='Create overlapping DTM tiles for mosaic')
parser.add_argument('--index_grid', metavar='path', required=True,
help='path to input overlapping index grid feature class')
parser.add_argument('--tiles_directory', metavar='path', required=True,
help='path to output DTM tiles directory')
args = parser.parse_args()
dtm_multi(index_grid=args.index_grid,
tiles_directory=args.tiles_directory)
... View more
02-04-2018
10:07 AM
|
0
|
1
|
886
|
POST
|
I'm able to create a reasonable index grid using the Grid Index Features tool based on the extent of my study area. What I'm trying to figure out is how to create an index grid that has a percentage of overlap. The following is going to be used to generate a DTM using Topo to Raster . The source is contours for the entire Dar es Salaam at a 1m resolution from surveyed contours. When I last tried to run it using the full contour dataset as a single feature class it ran for 7 days at a 5m resolution. Index Grid: Overlapping Index Grid: 10% DTM: 5m resolution
... View more
01-31-2018
08:30 AM
|
0
|
2
|
1461
|
POST
|
I have a road edge polygon that I'd like to split at each intersection. I came across the following post: How to buffer a polyline without overlap at verticies where Richard Fairhurst suggested a workflow to achieve the following. Unfortunately I wasn't able to follow his steps from where he created a buffer around the intersection points of the road centrelines. I've attached a sample of my road centrelines as well as my current buffer that I created based on the road centrelines and field that represents the road width. Any suggestions in how to achieve a split polygon based on the result that Richard Fairhurst achieved would be appreciated. from Richard Fairhurst reply Sample Dataset: (Road Centreline; Road Buffer (Merged); Road Intersection points)
... View more
12-24-2017
10:19 AM
|
0
|
0
|
811
|
POST
|
Thanks Dan , will check it out and get back to you.
... View more
12-19-2017
11:29 AM
|
0
|
0
|
2427
|
POST
|
I'm trying to determine how to calculate the Dice Similarity Coefficient between two rasters. The one raster is the ground truth result of a road surface area, the second raster is the result from a Computer Vision and Machine Learning (Convolutional Neural Network). I'm trying to understand how the Dice Similarity Coefficient works so that I can replicate it as a Python\NumPy function to run against all the sample results and ground truth results that were manually determined. NB: Id also like determine the following raster below which depicts True Positive; False Positive and False Negative Wikipedia: Sorensen-Dice Cofficient Formula: Raster Result:
... View more
12-19-2017
11:10 AM
|
0
|
2
|
10092
|
POST
|
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
09-23-2017
04:26 AM
|
1
|
3
|
3987
|
POST
|
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
09-19-2017
10:41 AM
|
0
|
5
|
4004
|
POST
|
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
09-19-2017
10:12 AM
|
0
|
0
|
4004
|
POST
|
I'm busy writing a Python module to calculate the NDVI (Normalized Difference Vegetation Index) based on the following post: Using Python to calculate NDVI with multiband imagery. According to the USGS background data on Landsat 8 Product, the tiles are provided as Digitial Numbers 16-bit unsigned integer format. The current examples (tutorials and reference material) that I've found to calculate NDVI use Landsat 5 + 7 which is stored as 8-bit unsigned integer format. Can I still use the same formula to calculate NDVI for Landsat 8 that is now stored as 16-bit unsigned integer format. The formula that I'm using is found in: Making Spatial Decisions Using GIS and Remote Sensing: A Workbook NDVI = (IR-R)/(IR+R) ArcGIS: USGS Landsat 8 GeoTIFF (Bands 4 + 5) ArcGIS: scale range of Landsat 8 Bands 4 + 5 USGS: Using the USGS Landsat 8 Product USGS EarthExplorer: Landsat 8 Dataset Selection USGS EarthExplorer: Landsat 8 Download Options
... View more
09-19-2017
09:41 AM
|
0
|
9
|
12217
|
POST
|
Hi Richard Thanks I clean forgot about Python's Pillow module. I was able to batch all the images and correct them.
... View more
09-18-2017
11:13 AM
|
0
|
0
|
825
|
POST
|
I have received imagery from a surveyor that is currently falling within the Northern Hemisphere and upside down. The following is easy to resolve if its a CAD drawing using the "Transformations" tab under the properties window of the CAD layer within ArcMap. Is there a way to apply the same transformation using Python for each image tile? CAD: Transformation Original Position: Northern Hemisphere and rotated. CAD: Transformation and Rotation CAD: Transformed (correct location) Imagery: Current Imagery (Northern hemisphere and rotated) Any assistance to achieve the following will be appreciated.
... View more
09-08-2017
01:43 AM
|
0
|
5
|
1565
|
Title | Kudos | Posted |
---|---|---|
3 | 01-16-2012 02:34 AM | |
1 | 05-07-2016 03:04 AM | |
1 | 04-10-2016 01:09 AM | |
1 | 03-13-2017 12:27 PM | |
1 | 02-17-2016 02:34 PM |
Online Status |
Offline
|
Date Last Visited |
03-04-2021
12:50 PM
|