|
POST
|
Hi Dan Thanks for your assistance once again, truly appreciated. I'm not sure Frequency Analysis is what I looking for. The COUNT_NAME field already has the number of buildings determined for each particular time interval. If you take "Asla Park B" as an example: Town Settlement Name TIME Field ( 15 - 25 min) TIME Field ( 20 - 25 min) Mossel Bay Asla Park B 6 Buildings = COUNT_NAME 148 Buildings = COUNT_NAME Let me kow if I've not explained it clearly enough. Thanks again
... View more
03-30-2016
02:21 PM
|
0
|
1
|
2198
|
|
POST
|
I have generated a summary statistics table based on service area analysis areas (5 min, 10 min, 15min etc) for each settlement area that I'm determining travel times from individual buildings within a settlement to a social facilities (i.e. Clinics, Schools etc) My output from my summary statistics table has the following columns: TOWN = Town Name SETTLEMENTNAME = Settlement Name NAME = Social Facility Name (i.e. School Name or Health Care Name) TIME = Service Area Time Interval ( 5 min, 10 min, 15min - 60 min) COUNT_NAME = The number of buildings within each service area (5 min, 10 min - 60 min) I'd like to convert the following Summary Statistics results table into the following format: TOWN SETTLEMENTNAME NAME 0 - 5 MIN = The number of buildings within 0 - 5 min of the Social Facility 5 - 10 MIN = The number of buildings within 5 - 10 min of the Social Facility 10 - 15 MIN = The number of buildings within 10 - 15 min of the Social Facility 15 - 20 MIN = The number of buildings within 15 - 20 min of the Social Facility 20 - 25 MIN = The number of buildings within 20 - 25 min of the Social Facility 25 - 30 MIN = The number of buildings within 25 - 30 min of the Social Facility 30 - 60 MIN = The number of buildings within 30- 60 min of the Social Facility I've attached a print screen of the Summary Statistics Table as well as the Output Results Table that I'm trying to achieve. I've briefly read that either Numpy or Pandas would be able to achieve the results that I'm looking for, but I'm a novice with Numpy and have never used Pandas. Any help would be appreciated.
... View more
03-30-2016
01:53 PM
|
0
|
3
|
6037
|
|
POST
|
Hi Dan I found that the new script that I wrote using Numpy has a further advantage in that it supports reclassifying an integer raster to float raster. '''
Created on Mar 19, 2016
@author: PeterW
'''
# import site-packages and modules
import os
import numpy as np
import arcpy
from arcpy.sa import *
# set arguments
input_raster = arcpy.Raster(r"E:\Projects\2016\G112030\SAA\Rasters\euc_road_clip")
output_folder = r"E:\Projects\2016\G112030\SAA\Rasters"
# set environment settings
arcpy.env.overwriteOutput = True
arcpy.env.outputCoordinateSystem = input_raster
# check out extensions
arcpy.CheckOutExtension("Spatial")
def reclass_raster(input_raster, output_folder):
lower_left = arcpy.Point(input_raster.extent.XMin,
input_raster.extent.YMin)
cell_size = input_raster.meanCellWidth
saa_array = arcpy.RasterToNumPyArray(input_raster)
reclass_array = np.zeros_like(saa_array)
bins = [0, 10, 25, 50, 100, 13069.5234375]
new_bins = [1, 1.25, 1.5, 2, 5]
new_classes = zip(bins[:-1], bins[1:], new_bins)
print(new_classes)
for rc in new_classes:
q1 = (saa_array >= rc[0])
q2 = (saa_array < rc[1])
z = np.where(q1 & q2, rc[2], 0)
reclass_array = reclass_array + z
reclass_raster = arcpy.NumPyArrayToRaster(reclass_array,
lower_left,
x_cell_size=cell_size,
value_to_nodata=0)
reclass_raster.save(os.path.join(output_folder, "cost_road"))
reclass_raster(input_raster, output_folder)
# check in extensions
arcpy.CheckInExtension("Spatial")
... View more
03-30-2016
03:23 AM
|
0
|
1
|
2782
|
|
POST
|
Hi Steve I was originally using map algebra Conditional Statement to reclassify my output raster and was looking for a way to build a function for the new classes that I wanted the raster to be reclassified to. Darren and Dan suggested that Numpy would allow me to generate a function to define the classes more cleanly. I must admit I had not considered the reclassify tool, RemapRange would have worked just just as well.
... View more
03-29-2016
01:24 AM
|
0
|
6
|
2782
|
|
POST
|
Hi Dan Your sample code worked perfectly, with some minor changes. '''
Created on Mar 19, 2016
@author: PeterW
'''
# import site-packages and modules
import os
import numpy as np
import arcpy
from arcpy.sa import *
# set arguments
input_raster = arcpy.Raster(r"E:\Python\Testing\Numpy_Reclassify\SAA_Raster\tbr")
output_folder = r"E:\Python\Testing\Numpy_Reclassify\SAA_Reclassified"
# set environment settings
arcpy.env.overwriteOutput = True
arcpy.env.outputCoordinateSystem = input_raster
# check out extensions
arcpy.CheckOutExtension("Spatial")
def reclass_raster(input_raster, output_folder):
lower_left = arcpy.Point(input_raster.extent.XMin,
input_raster.extent.YMin)
cell_size = input_raster.meanCellWidth
saa_array = arcpy.RasterToNumPyArray(input_raster)
reclass_array = np.zeros_like(saa_array)
bins = [0, 5, 10, 15, 20, 25, 30, 60]
new_bins = [5, 10, 15, 20, 25, 30, 60]
new_classes = zip(bins[:-1], bins[1:], new_bins)
for rc in new_classes:
q1 = (saa_array >= rc[0])
q2 = (saa_array < rc[1])
z = np.where(q1 & q2, rc[2], 0)
reclass_array = reclass_array + z
reclass_raster = arcpy.NumPyArrayToRaster(reclass_array,
lower_left,
x_cell_size=cell_size,
value_to_nodata=0)
Int(reclass_raster).save(os.path.join(output_folder, "rec"))
reclass_raster(input_raster, output_folder)
# check in extensions
arcpy.CheckInExtension("Spatial") Regards
... View more
03-28-2016
02:21 PM
|
0
|
1
|
3387
|
|
POST
|
Hi Dan The following looks to be what I was looking for. I've been swamped the last few days. Will come back to the following over the weekend and post my code using numpy.
... View more
03-17-2016
01:11 AM
|
1
|
0
|
3387
|
|
POST
|
I've written a Python script that reclassifies a raster using a set of conditional statements into 4 classes: 0 - 12 12 - 24 24 - 30 30 - 60 rec_ras = Con(costmin <= 12, 12, Con((costmin > 12) & (costmin <= 24), 24, Con((costmin > 24) & (costmin <= 30), 30, Con((costmin > 30) & (costmin <= 60), 60)))) I then needed to change the classes into 7 classes: 0 - 5 5 - 10 10 - 15 15 - 20 20 - 25 25 - 30 30 - 60 rec_ras = Con(costmin <= 5, 5, Con((costmin > 5) & (costmin <= 10), 10, Con((costmin > 10) & (costmin <= 15), 15, Con((costmin > 15) & (costmin <= 20), 20, Con((costmin > 20) & (costmin <= 25), 25, Con((costmin > 25) & (costmin <= 30), 30, Con((costmin > 30) & (costmin <= 60), 60))))))) Is ther a better way of creating a Python function that can create my conditional statement as and if when I need to change the classes before running it. The following is part of a larger Python Class.
... View more
03-11-2016
10:42 AM
|
0
|
20
|
9995
|
|
POST
|
Hi Jason After fixing the line 98 incorrect allocation of fac (Flow Accumulation) to Flow Direction, my script completes within the same amount of time as my previous script. I'm going to look into developing a function to allocate the variables automatically as this is standardised. The only input required from the user is the: Input DEM Raster Folder File Geodatabase
... View more
02-27-2016
10:05 AM
|
0
|
1
|
2402
|
|
POST
|
Hi Jason Thanks for your feedback, its truly appreciated as I'm still a novice when it comes to programming. With regards to your suggestions. Would you mind giving me an example to explain your first suggestion: Group all of your methods at the top under imports instead of co-mingling functions and statements. I only realized once you pointed it out that I was defining my variables multiple times, thanks for pointing it out. Lastly could you give me some advice how I could create a function that would build my input and output arguments based on the raster workspace and fgdb.
... View more
02-26-2016
12:29 PM
|
0
|
1
|
2402
|
|
POST
|
Hi Luke If you look at my original Python Script I used os.path.join and decided to move to PathLib Path. The unfortunate thing is that all ArcPy scripts require paths to be string and not objects hence using strings to create the paths from the Path objects.
... View more
02-25-2016
10:06 AM
|
1
|
5
|
3912
|
|
POST
|
Hi Dan Thanks for the following, I'll have a look at the following this evening and get back to you if I get stuck, if that fine with you.
... View more
02-25-2016
05:32 AM
|
1
|
7
|
3912
|
|
POST
|
I'm busy rewriting my hydrological model python scripts into functions that are modular and found that my original python script runs twice as fast as my updated python script. The only reason I can find is that the original scripts outputs are saved in memory before being passed onto the next function. My new python script is writing and reading each output from disk. I might be incorrect in my understanding. If anyone can give me guidance on improving my new python script in both pythonic structure and performance. Original Python Script '''
Created on May 20, 2015
@author: PeterW
'''
# import system modules and site packages
import os
import time
import arcpy
import ArcHydroTools
# check out Spatial Analyst Extension
arcpy.CheckOutExtension("Spatial")
# set environment settings
arcpy.env.overwriteOutput = True
# set input and output arguments
raw = r"E:\Python\Temp\Model04\DEM\raw"
rasWs = r"E:\Python\Temp\Model04\Layers04"
outWs = r"E:\Python\Temp\Model04\Model04.gdb"
# Processing time
def hms_string(sec_elapsed):
h = int(sec_elapsed / (60 * 60))
m = int(sec_elapsed % (60 * 60) / 60)
s = sec_elapsed % 60
return "{}h:{:>02}m:{:>05.2f}s".format(h, m, s)
start_time1 = time.time()
# archydro variables
fill_sinks = os.path.join(rasWs, "fil")
flow_dir = os.path.join(rasWs, "fdr")
flow_acc = os.path.join(rasWs, "fac")
streams = os.path.join(rasWs, "str")
stream_seg = os.path.join(rasWs, "strlnk")
catchment_grid = os.path.join(rasWs, "cat")
catchment_poly = os.path.join(outWs, "Layers","Catchment")
drainage_line = os.path.join(outWs, "Layers", "DrainageLine")
adj_catch = os.path.join(outWs, "Layers", "AdjointCatchment")
try:
# calculate the fill sinks
arcpy.AddMessage("Processing Fill Sinks")
ArcHydroTools.FillSinks(raw, fill_sinks)
# calculate the flow direction
arcpy.AddMessage("Processing Flow Direction")
ArcHydroTools.FlowDirection(fill_sinks, flow_dir)
# calculate the flow accumulation
arcpy.AddMessage("Processing Flow Accumulation")
ArcHydroTools.FlowAccumulation(flow_dir, flow_acc)
# calculate the maximum flow accumulation
arcpy.AddMessage("Processing Flow Accumulation Maximum")
maxcellsResult = arcpy.GetRasterProperties_management(flow_acc, "MAXIMUM")
maxcells = maxcellsResult.getOutput(0)
print maxcells
# calculate the stream threshold number of cells
arcpy.AddMessage("Processing Stream Threshold")
stream_threshold_numcells = (int(maxcells)*0.125/100)
print stream_threshold_numcells
# calculate the stream definition
arcpy.AddMessage("Processing Stream Definition")
ArcHydroTools.StreamDefinition(flow_acc, stream_threshold_numcells, streams)
# calculate the stream segmentation
arcpy.AddMessage("Processing Stream Segmentation")
ArcHydroTools.StreamSegmentation(streams, flow_dir, stream_seg)
# calculate the catchment grid delineation
arcpy.AddMessage("Processing Catchment Grid Delineation")
ArcHydroTools.CatchmentGridDelineation(flow_dir, stream_seg, catchment_grid)
# calculate the catchment polygons from the catchment grid
arcpy.AddMessage("Processing Catchment Polygons")
ArcHydroTools.CatchmentPolyProcessing(catchment_grid, catchment_poly)
# calculate the drainage lines from the stream segmentation grid
arcpy.AddMessage("Processing DrainageLines")
ArcHydroTools.DrainageLineProcessing(stream_seg, flow_dir, drainage_line)
# calculate the adjoint catchment polygons
arcpy.AddMessage("Processing Ajdoint Catchments")
ArcHydroTools.AdjointCatchment(drainage_line, catchment_poly, adj_catch)
arcpy.AddMessage("Completed Processing archydro Main Model")
except:
print(arcpy.GetMessages(2))
pass
# Determine the time take to process hydrological characteristics
end_time1 = time.time()
print ("It took {} to process hydrological characteristics".format(hms_string(end_time1 - start_time1)))
arcpy.CheckInExtension("Spatial") New Python Script '''
Created on Feb 24, 2016
@author: PeterW
'''
# import system modules and site packages
import time
from pathlib import Path
import arcpy
import ArcHydroTools
# check out extension
arcpy.CheckOutExtension("Spatial")
# set environment settings
arcpy.env.overwriteOutput = True
# set input and output arguments
dem = r"E:\Python\Temp\Model04\DEM\raw"
raster_workspace = r"E:\Python\Temp\Model04\Layers04"
fgdb = Path(r"E:\Python\Temp\Model04\Model04.gdb")
# Processing time
def hms_string(sec_elapsed):
h = int(sec_elapsed / (60 * 60))
m = int(sec_elapsed % (60 * 60) / 60)
s = sec_elapsed % 60
return "{}h:{:>02}m:{:>05.2f}s".format(h, m, s)
start_time1 = time.time()
# generate fill sinks grid
def fill_sinks(dem, raster_workspace):
fill_sinks = "{0}\\{1}".format(raster_workspace, "fil")
arcpy.AddMessage("Processing Fill Sinks")
ArcHydroTools.FillSinks(dem, fill_sinks)
fill_sinks(dem, raster_workspace)
# generate flow direction grid
def flow_direction(raster_workspace):
fill_sinks = "{0}\\{1}".format(raster_workspace, "fil")
flow_direction = "{0}\\{1}".format(raster_workspace, "fdr")
arcpy.AddMessage("Processing Flow Direction")
ArcHydroTools.FlowDirection(fill_sinks, flow_direction)
flow_direction(raster_workspace)
# generate flow accumulation grid
def flow_accumulation(raster_workspace):
flow_direction = "{0}\\{1}".format(raster_workspace, "fdr")
flow_accumulation = "{0}\\{1}".format(raster_workspace, "fac")
arcpy.AddMessage("Processing Flow Accumulation")
ArcHydroTools.FlowAccumulation(flow_direction, flow_accumulation)
flow_accumulation(raster_workspace)
# calculate stream threshold based on 0.5% of maximum flow accumulation
def stream_threshold(raster_workspace):
flow_accumulation = "{0}\\{1}".format(raster_workspace, "fac")
arcpy.AddMessage("Processing Flow Stream Threshold")
maxcellsResult = arcpy.GetRasterProperties_management(flow_accumulation, "MAXIMUM")
maxcells = maxcellsResult.getOutput(0)
arcpy.AddMessage("{} Maximum Cells".format(maxcells))
stream_threshold_numcells = (int(maxcells)*0.5/100)
arcpy.AddMessage("{} Stream Threshold".format(stream_threshold_numcells))
return stream_threshold_numcells
numcells = stream_threshold(raster_workspace)
# generate the stream definition grid
def stream_definition(raster_workspace):
flow_accumulation = "{0}\\{1}".format(raster_workspace, "fac")
stream = "{0}\\{1}".format(raster_workspace, "str")
arcpy.AddMessage("Processing Stream Definition")
ArcHydroTools.StreamDefinition(flow_accumulation, numcells, stream)
stream_definition(raster_workspace)
# generate the stream segmentation grid
def stream_segmentation(raster_workspace):
stream = "{0}\\{1}".format(raster_workspace, "str")
flow_direction = "{0}\\{1}".format(raster_workspace, "fac")
stream_link = "{0}\\{1}".format(raster_workspace, "strlnk")
arcpy.AddMessage("Processing Stream Segmentation")
ArcHydroTools.StreamSegmentation(stream, flow_direction, stream_link)
stream_segmentation(raster_workspace)
# calculate the catchment grid delineation
def catchment_grid(raster_workspace):
flow_direction = "{0}\\{1}".format(raster_workspace, "fdr")
stream_link = "{0}\\{1}".format(raster_workspace, "strlnk")
catchment_grid = "{0}\\{1}".format(raster_workspace, "cat")
arcpy.AddMessage("Processing Catchment Grid Delineation")
ArcHydroTools.CatchmentGridDelineation(flow_direction, stream_link, catchment_grid)
catchment_grid(raster_workspace)
# calculate the catchment polygons from the catchment grid
def catchment_polygon(raster_workspace, fgdb):
catchment_grid = "{0}\\{1}".format(raster_workspace, "cat")
catchment_polygon = "{0}\\{1}".format(Path(fgdb, "Layers"), "Catchment")
arcpy.AddMessage("Processing Catchment Polygon")
ArcHydroTools.CatchmentPolyProcessing(catchment_grid, catchment_polygon)
catchment_polygon(raster_workspace, fgdb)
# calculate the drainage lines from the stream segmentation grid
def drainage_line(raster_workspace, fgdb):
stream_link = "{0}\\{1}".format(raster_workspace, "strlnk")
flow_direction = "{0}\\{1}".format(raster_workspace, "fdr")
drainage_line = "{0}\\{1}".format(Path(fgdb, "Layers"), "DrainageLine")
arcpy.AddMessage("Processing DrainageLine")
ArcHydroTools.DrainageLineProcessing(stream_link, flow_direction, drainage_line)
drainage_line(raster_workspace, fgdb)
# calculate the adjoint catchment polygons
def adjoint_catchment(fgdb):
drainage_line = "{0}\\{1}".format(Path(fgdb, "Layers"), "DrainageLine")
catchment_polygon = "{0}\\{1}".format(Path(fgdb, "Layers"), "Catchment")
adjoint_catchment = "{0}\\{1}".format(Path(fgdb, "Layers"), "AdjointCatchment")
arcpy.AddMessage("Processing Adjoint Catchment")
ArcHydroTools.AdjointCatchment(drainage_line, catchment_polygon, adjoint_catchment)
adjoint_catchment(fgdb)
arcpy.CheckInExtension("Spatial")
# Determine the time take to process hydrological characteristics
end_time1 = time.time()
print ("It took {} to process hydrological characteristics".format(hms_string(end_time1 - start_time1)))
... View more
02-24-2016
03:58 PM
|
0
|
16
|
10356
|
|
POST
|
Hi Paul I've not heard back from you, have you come right with the following?
... View more
02-24-2016
02:18 AM
|
0
|
0
|
647
|
|
POST
|
Hi Alec If you can send me an LinkedIn invite @ Peter Wilson I'll accept and send you my details. I don't want to post my work email address on the net. Regards
... View more
02-21-2016
02:30 PM
|
0
|
0
|
647
|
|
POST
|
Hi Alec I'm a GIS Engineer at Aurecon South Africa, Cape Town and I specialize in hydrological and hydraulic modelling. If you could send me shapefile or kml file of the extent of your study area, I'd like to see if I can help you out.
... View more
02-21-2016
02:33 AM
|
1
|
3
|
3366
|
| 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
|