/
raster_function.py
133 lines (119 loc) · 5.62 KB
/
raster_function.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
import arcpy
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")
import numpy as np
import pandas as pd
import gc
#####################################################################################################################
def catcsv2raster(inCSV, Value, inTemplate, outRaster, dtype='Int', idName='COMID'):
'''
__author__ = "Ryan Hill <hill.ryan@epa.gov>"
"Marc Weber <weber.marc@epa.gov>"
Converts CSV table to GeoTIFF and save output.
Requires template raster with NHDPlusV2 COMIDs that match CSV files.
Arguments
---------
inCSV : csv table
Input catchment values
Value : Text input from user
Column in table to use as values
inTemplate : GeoTIFF
Template raster with COMIDs for pixel values
outRaster : GeoTIFF
Output raster file path and name
dtype : Text input from user
Desired data type of output raster - default = None
If no dtype is provided, the default is to produce an 'Int' raster
idName : Text input from user
Column in table to use as unique ID - default = 'COMID'
'''
#Read raster and CSV
inRas = arcpy.Raster(inTemplate)
#Fix environment settings to inRas
dsc=arcpy.Describe(inRas)
arcpy.env.extent=dsc.Extent
arcpy.env.compression = 'LZW'
ext=dsc.Extent
ll = arcpy.Point(ext.XMin, ext.YMin)
arcpy.env.outputCoordinateSystem=dsc.SpatialReference
arcpy.env.cellSize=dsc.meanCellWidth
arcpy.env.mask = inRas
cellSize = dsc.meanCellWidth
rstArray = arcpy.RasterToNumPyArray(inRas)
lookup = pd.read_csv(inCSV)
#Prep data by adding noData number to front of vectors
b = np.float64(np.append(-9999, np.array(lookup[idName])))
c = np.float64(np.append(-9999, np.array(lookup[Value])))
a = rstArray.flatten()
a[a==0] = -9999
a = np.where(np.in1d(a,b), a, -9999)
#Run numpy query to replace COMID raster with desired values:
bsort = np.argsort(b) #Create sorting index
apos = np.searchsorted(b[bsort], a) #Search a on sorted b
indices = bsort[apos] #Get indices in b that match a
z = c[indices] #Make final vector from desired data (c)
z.shape = rstArray.shape #Reshape back to 2d
newRas = arcpy.NumPyArrayToRaster(z, lower_left_corner=ll, x_cell_size=cellSize, y_cell_size=cellSize, value_to_nodata=-9999)
if dtype == 'Int':
arcpy.CopyRaster_management(newRas, outRaster, "", "", "", "", "", "16_BIT_SIGNED")
else:
arcpy.CopyRaster_management(newRas, outRaster, "", "", "", "", "", "32_BIT_FLOAT")
del newRas, a, b, c, bsort, apos, indices, z, inRas, rstArray
gc.collect()
#####################################################################################################################
#####################################################################################################################
def catcsv2raster2(lookup, Value, inTemplate, outRaster, dtype='Int', idName='COMID'):
'''
__author__ = "Ryan Hill <hill.ryan@epa.gov>"
"Marc Weber <weber.marc@epa.gov>"
Converts CSV table to GeoTIFF and save output.
Requires template raster with NHDPlusV2 COMIDs that match CSV files.
Arguments
---------
lookup : pandas table
Input catchment values
Value : Text input from user
Column in table to use as values
inTemplate : GeoTIFF
Template raster with COMIDs for pixel values
outRaster : GeoTIFF
Output raster file path and name
dtype : Text input from user
Desired data type of output raster - default = None
If no dtype is provided, the default is to produce an 'Int' raster
idName : Text input from user
Column in table to use as unique ID - default = 'COMID'
'''
#Read raster and CSV
inRas = arcpy.Raster(inTemplate)
#Fix environment settings to inRas
dsc=arcpy.Describe(inRas)
arcpy.env.extent=dsc.Extent
arcpy.env.compression = 'LZW'
ext=dsc.Extent
ll = arcpy.Point(ext.XMin, ext.YMin)
arcpy.env.outputCoordinateSystem=dsc.SpatialReference
arcpy.env.cellSize=dsc.meanCellWidth
arcpy.env.mask = inRas
cellSize = dsc.meanCellWidth
rstArray = arcpy.RasterToNumPyArray(inRas)
# Prep data by adding noData number to front of vectors
b = np.float64(np.append(-9999, np.array(lookup[idName])))
c = np.float64(np.append(-9999, np.array(lookup[Value])))
a = rstArray.flatten()
a[a==0] = -9999
a = np.where(np.in1d(a,b), a, -9999)
#Run numpy query to replace COMID raster with desired values:
bsort = np.argsort(b) #Create sorting index
apos = np.searchsorted(b[bsort], a) #Search a on sorted b
indices = bsort[apos] #Get indices in b that match a
z = c[indices] #Make final vector from desired data (c)
z.shape = rstArray.shape #Reshape back to 2d
newRas = arcpy.NumPyArrayToRaster(z, lower_left_corner=ll, x_cell_size=cellSize, y_cell_size=cellSize, value_to_nodata=-9999)
if dtype == 'Int':
arcpy.CopyRaster_management(newRas, outRaster, "", "", "", "", "", "16_BIT_SIGNED")
else:
arcpy.CopyRaster_management(newRas, outRaster, "", "", "", "", "", "32_BIT_FLOAT")
del newRas, a, b, c, bsort, apos, indices, z, inRas, rstArray
gc.collect()
#####################################################################################################################