· 8 years ago · Apr 11, 2018, 11:14 AM
1python format_veg_params.py <template raster> <land cover raster>
2<output veg parameter file> <land cover classification scheme>
3
4Traceback (most recent call last):
5 File "format_veg_params.py", line 180, in <module>
6 main()
7 File "format_veg_params.py", line 174, in main
8 format_veg_params(sys.argv[1],sys.argv[2],sys.argv[3],sys.argv[4])
9 File "format_veg_params.py", line 120, in format_veg_params
10 clsAttributes)-1)].ravel()).argmax()
11TypeError: Cannot cast array data from dtype('float32') to dtype('int64') according to the rule 'safe'
12
13# import dependencies
14import os
15import sys
16import json
17import numpy as np
18from osgeo import gdal
19from osgeo.gdalnumeric import *
20from osgeo.gdalconst import *
21
22def format_veg_params(basinMask,lcData,outVeg,scheme='IGBP'):
23 """
24 FUNCTION: format_veg_params
25 ARGUMENTS: basinMask - path to basin template raster
26 lcData - path to land cover raster
27 outveg - path output vegetation parameter file
28 KEYWORDS: scheme - Abbreviation of land cover classification scheme the
29 input land cover data is formatted in
30 RETURNS: n/a
31 NOTES: Returns no variables but writes an output file
32 """
33
34 # define script file path for relative path definitions
35 __location__ = os.path.realpath(
36 os.path.join(os.getcwd(), os.path.dirname(__file__)))
37
38 # get land cover classification scheme and create path to lookup table
39 if scheme == 'IGBP':
40 attriFile = os.path.join(__location__,'veg_type_attributes_igbp.json')
41 elif scheme == 'GLCC':
42 attriFile = os.path.join(__location__,'veg_type_attributes_glcc.json')
43 elif scheme == 'IPCC':
44 attriFile = os.path.join(__location__,'veg_type_attributes_ipcc.json')
45 else:
46 raise SyntaxError('Land cover classification scheme not supported')
47
48 band = 1 # constant variable for reading in data
49
50 # open/read veg scheme json file
51 with open(attriFile) as data_file:
52 attriData = json.load(data_file)
53
54 # pass look up information into variable
55 clsAttributes = attriData['classAttributes']
56
57 # create list of input raster files
58 infiles = [os.path.join(__location__,basinMask),
59 os.path.join(__location__,lcData)]
60
61 try: # try to read in the raster data
62
63 # read basin grid raster
64 ds = gdal.Open(infiles[0],GA_ReadOnly)
65 b1 = ds.GetRasterBand(band)
66 mask = BandReadAsArray(b1)
67 maskRes = ds.GetGeoTransform()[1]
68 ds = None
69 b1 = None
70
71 # read land cover raster
72 ds = gdal.Open(infiles[1],GA_ReadOnly)
73 b1 = ds.GetRasterBand(band)
74 lccls = BandReadAsArray(b1)
75 clsRes = ds.GetGeoTransform()[1]
76 ds = None
77 b1 = None
78
79 # if not working, give error message
80 except AttributeError:
81 raise IOError('Raster file input error, check that all paths are correct')
82
83 ratio = maskRes/clsRes # get ratio of high resoltion to low resolution
84
85 # get file path to output file
86 vegfile = os.path.join(__location__,outVeg)
87
88 # check if the output parameter file exists, if so delete it
89 if os.path.exists(vegfile)==True:
90 os.remove(vegfile)
91
92 try: # try to write output veg parameter file
93
94 # open output file for writing
95 with open(vegfile, 'w') as f:
96
97 cnt = 1 # grid cell id counter
98
99 # loop over each pixel in the template raster
100 for i in range(mask.shape[0]):
101 y1 = int(i*ratio)
102 y2 = int(y1+ratio)
103 for j in range(mask.shape[1]):
104 x1 = int(j*ratio)
105 x2 = int(x1+ratio)
106
107 # get land cover data within template raster pixel
108 tmp = lccls[y1:y2,x1:x2]
109
110 # if not a mask value...
111 if mask[i,j] == 1:
112 # if there are nodata values...
113 if np.any(tmp>len(clsAttributes)-1)==True:
114 # ...find where they are...
115 negdx = np.where(tmp>len(clsAttributes)-1)
116 try:
117 # ...and calculate hisogram for only data values
118 tmp[negdx] = np.bincount(tmp[np.where(tmp<len(
119 clsAttributes)-1)].ravel()).argmax()
120 except ValueError:
121 tmp[negdx] = 0
122 uniqcnt = np.unique(tmp) # number of unique values
123 clscnt = np.bincount(tmp.ravel())
124
125 # check if there is only 1 value in unique value array
126 if type(uniqcnt).__name__ == 'int':
127 Nveg = 1 # if so, only 1 value
128 else:
129 Nveg = len(uniqcnt) # if not, n veg equal to unique values
130
131 # write the grid cell id and number of classes
132 f.write('{0} {1}n'.format(cnt,Nveg))
133
134 # if there is vegetation data...
135 if Nveg != 0:
136 # ...loop over each class...
137 for t in range(uniqcnt.size):
138 # ...and grab the parameter attributes
139 vegcls = int(uniqcnt[t]) # class value
140 Cv = np.float(clscnt[uniqcnt[t]])/np.float(clscnt.sum()) # percent coverage
141 attributes = clsAttributes[vegcls]['properties'] # intermediate variale
142 rdepth1 = str(attributes['rootd1']) # rooting depth layer 1
143 rfrac1 = str(attributes['rootfr1']) # rooting fraction layer 1
144 rdepth2 = str(attributes['rootd2']) # rooting depth layer 2
145 rfrac2 = str(attributes['rootfr2']) # rooting fraction layer 2
146 rdepth3 = str(attributes['rootd3']) # rooting depth layer 3
147 rfrac3 = str(attributes['rootfr3']) # rooting fraction layer 3
148
149 # write the veg class information
150 f.write('t{0} {1:.4f} {2} {3} {4} {5} {6} {7}n'.format(vegcls,Cv,rdepth1,rfrac1,rdepth2,rfrac2,rdepth3,rfrac3))
151
152 cnt+=1 # plus one to grid cell id counter
153
154 # except raise an error when it doesn't work
155 except IOError:
156 raise IOError('Cannot write output file, error with output veg parameter file path')
157
158 return
159
160def main():
161 n_args = len(sys.argv)
162
163 # Check user inputs
164 if n_args != 5:
165 print("Wrong user input")
166 print("Script writes the vegetation parameter file for the VIC model")
167 print("usage: python format_veg_params.py <template raster> <land cover raster> <output veg parameter file> <land cover classification scheme>")
168 print("Exiting system...")
169 sys.exit()
170
171 else:
172 # Pass command line arguments into function
173 format_veg_params(sys.argv[1],sys.argv[2],sys.argv[3],sys.argv[4])
174
175 return
176
177# Execute the main level program if run as standalone
178if __name__ == "__main__":
179 main()