-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathburn_raster_driver_b.py
279 lines (213 loc) · 13.3 KB
/
burn_raster_driver_b.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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
#burn raster driver
import os
import csv
import pandas as pd
import gdal
import osr
import subprocess as sub
import sys
import shutil
import utm
#plan: use processed csv to rasterize a new source raster, and run the raster burner hack of the chi_mapping_tool
#directory = 'R:\\LSDTopoTools\\Topographic_projects\\full_himalaya\\'
directory = '/exports/csce/datastore/geos/users/s1134744/LSDTopoTools/Topographic_projects/full_himalaya/'
#data_source = '/exports/csce/datastore/geos/users/s1134744/LSDTopoTools/Topographic_projects/TRMM_data/full_data/non_monsoon_trans.tif'
#data_source = 'R:\\exhumation\\0_2.tif'#source of the rasterized exhumation data
data_source_exhumation = '/exports/csce/datastore/geos/users/s1134744/exhumation/0_2.tif'
data_source_glaciated = '/exports/csce/datastore/geos/users/s1134744/LSDTopoTools/Topographic_projects/shapefiles/glims_ice/ice_100.bil'
data_source_distance = '/exports/csce/datastore/geos/users/s1134744/LSDTopoTools/Topographic_projects/shapefiles/simplified_distance_metres_1km.tif'
data_source_distance_along = '/exports/csce/datastore/geos/users/s1134744/LSDTopoTools/Topographic_projects/shapefiles/euc_allocation_5.tif'
data_source_monsoon = '/exports/csce/datastore/geos/users/s1134744/LSDTopoTools/Topographic_projects/TRMM_data/jun_jul_aug.tif'
data_source_strain = '/exports/csce/datastore/geos/users/s1134744/LSDTopoTools/Topographic_projects/strain/second_invariant.tif'
data_source_tectonics = '/exports/csce/datastore/geos/users/s1134744/LSDTopoTools/Topographic_projects/shapefiles/fault_zones/digitized.bil'#source of digitized shapefile
def utmExtents(summary_directory,fname,corner):
with open (summary_directory+fname+corner+".txt", "r") as text:
data = text.read().replace(',', '')
data = data.replace('(','')
data = data.replace(')','')
data = data.split()
coordinates = [data[2],data[3]]
return coordinates
def utmToLatLon(coordinate,proj4):
zone = proj4.split(' ')
zone = zone[1]
zone = zone.replace('+zone=','')
easting = float(coordinate[0])
northing = float(coordinate[1])
lat_lon_coordinate = utm.to_latlon(int(easting),int(northing),int(zone),northern = True)
return lat_lon_coordinate
def getGlaciers(full_directory,summary_directory,fname,write_name,SRTM90=False):
#these are in UTM, needs to be lat lon
LL = utmExtents(summary_directory,fname,corner="_lower_left")
UR = utmExtents(summary_directory,fname,corner="_upper_right")
simplified_geology = '/exports/csce/datastore/geos/users/s1134744/LSDTopoTools/Topographic_projects/shapefiles/fault_zones/digitized.shp' #source of digitized shapefile
glaciers = '/exports/csce/datastore/geos/users/s1134744/LSDTopoTools/Topographic_projects/shapefiles/glims_ice/clipped_exists.shp'
ds = gdal.Open(summary_directory+fname+'.bil')
#getting target srs from current DEM
ds = ds.GetProjection()
ds = osr.SpatialReference(ds)
ds = ds.ExportToProj4()
try:
LL_lat_lon = utmToLatLon(LL,ds)
UR_lat_lon = utmToLatLon(UR,ds)
except:
sys.exit()
#clipping
try:
print('ogr2ogr -f "ESRI Shapefile" %s.shp %s -clipsrc %s %s %s %s' %(summary_directory+fname+'_glaciers',glaciers,LL_lat_lon[0],LL_lat_lon[1],UR_lat_lon[0],UR_lat_lon[1]))
sys.exit()
except:
sys.exit()
os.system(('ogr2ogr -f "ESRI Shapefile" %s.shp %s -clipsrc %s %s %s %s') %(summary_directory+fname+'_glaciers',glaciers,LL_lat_lon[0],LL_lat_lon[1],UR_lat_lon[0],UR_lat_lon[1]))
#transform crs and rasterize
os.system('ogr2ogr -t_srs'+" '"+ds+"' "+summary_directory+fname+'_glaciers_utm'+'.shp '+summary_directory+fname+'_glaciers.shp')
if os.path.isfile(full_directory+fname+'_glaciated.bil'):
try:
os.remove(full_directory+fname+'_glaciated.bil')
os.remove(full_directory+fname+'_glaciated.hdr')
except:
print("Error removing raster files. Forcing exit to prevent unseen errors")
sys.exit()
if SRTM90:
os.system('gdal_rasterize -of ENVI -a exists -a_nodata 0 -tr 90 90 -l '+fname+'_glaciers_utm'+' '+summary_directory+fname+'_glaciers_utm'+'.shp '+full_directory+write_name+'_glaciated'+'.bil')
else:
os.system('gdal_rasterize -of ENVI -a exists -a_nodata 0 -tr 30 30 -l '+fname+'_glaciers_utm'+' '+summary_directory+fname+'_glaciers_utm'+'.shp '+full_directory+write_name+'_glaciated'+'.bil')
def getRaster(full_directory,summary_directory,fname,write_name,raster_source,raster_name='to_burn',SRTM90=False):
LL = utmExtents(summary_directory,fname,corner="_lower_left")
UR = utmExtents(summary_directory,fname,corner="_upper_right")
#getting target srs from current DEM
try:
print full_directory+fname+'.bil'
ds = gdal.Open(full_directory+fname+'.bil')
ds = ds.GetProjection()
ds = osr.SpatialReference(ds)
ds = ds.ExportToProj4()
except:
print("error getting burn raster, skipping...")
return
#clipping raster, reprojecting, and setting resolution
res_90 = "gdalwarp -ot Float64 -of ENVI -tr 90 90 -s_srs 'EPSG:4326' -t_srs"+" '"+ds+"' -te %s %s %s %s %s %s%s_%s.bil" %(LL[0],LL[1],UR[0],UR[1],raster_source,full_directory,fname,raster_name)
res_30 = "gdalwarp -ot Float64 -of ENVI -tr 30 30 -s_srs 'EPSG:4326' -t_srs"+" '"+ds+"' -te %s %s %s %s %s %s%s_%s.bil" %(LL[0],LL[1],UR[0],UR[1],raster_source,full_directory,fname,raster_name)
#before rasterising, check if raster is already present in directory, remove if it is
if os.path.isfile(full_directory+fname+'_'+raster_name+'.bil'):
try:
os.remove(full_directory+fname+'_'+raster_name+'.bil')
os.remove(full_directory+fname+'_'+raster_name+'.hdr')
except:
print("Error removing raster files. Forcing exit to prevent unseen errors")
sys.exit()
print 'srtm 90 is ...s...', SRTM90
if SRTM90:
os.system(res_90)
print res_90
else:
os.system(res_30)
print res_30
#opening processed source file to access directory structure
name_list = ['himalaya_processed','himalaya_b_processed','himalaya_c_processed']
for name in name_list:
with open(directory+name+'.csv','r') as csvfile:
csvReader = csv.reader(csvfile,delimiter=',')
next(csvReader)
for row in csvReader:
#print(row)
#generating target path
#for writing output, will make it easier in handling the generated files with existing scripts
part_1 = str(row[0])
part_1 = part_1.replace('.','_')
part_2 = str(("%.2f" %float(row[2])))+'_'+str(("%.2f" %float(row[3])))
part_2 = part_2.replace('.','_')
full_target = os.path.join(directory+part_1,part_2+'_'+part_1+'_'+str(row[1]),str(row[5])+'/')
#print(full_target)
#summary directory. This is where the extent.txt files are preserved from the original analysis
summary_target = directory+part_1+'_summary/'
fname = part_1+'_'+str(row[1])
max_basin = (int(row[6])/2)+int(row[5])
write_name = str(row[1])+str(row[5])+'_'+str(max_basin)
#column_header = "digitized"
secondary_column_header = "Null"
if not os.path.isfile(full_target+fname+'.bil'):
try:
shutil.copy2(summary_target+fname+'.bil',full_target+fname+'.bil')
shutil.copy2(summary_target+fname+'.hdr',full_target+fname+'.hdr')
except:
print("cannot find the dem for %s"%(fname))
#adaptions for concavity catcher. Needs to get precipitation and litho rasters from summary directory
#shutil.copy2(summary_target+fname+'_precipitation.bil',full_target+write_name+'_precipitation.bil')
#shutil.copy2(summary_target+fname+'_precipitation.hdr',full_target+write_name+'_precipitation.hdr')
#shutil.copy2(summary_target+fname+'_geology.bil',full_target+write_name+'_geology.bil')
#shutil.copy2(summary_target+fname+'_geology.hdr',full_target+write_name+'_geology.hdr')
#^^^ not necessary as rasters are sourced from the summary directory
#getRaster(full_target,summary_target,fname,write_name,data_source_exhumation,raster_name='exhumation')
#getRaster(full_target,summary_target,fname,write_name,data_source_glaciated,raster_name='glaciated')
#getRaster(full_target,summary_target,fname,write_name,data_source_distance,raster_name='distance')
#getRaster(full_target,summary_target,fname,write_name,data_source_monsoon,raster_name='monsoon')
getRaster(full_target,summary_target,fname,write_name,data_source_strain,raster_name='strain_ezz')
#getGlaciers(full_target,summary_target,fname,write_name)
#print(full_target,summary_target,fname,write_name)
#try:
#getGeologyRaster(full_target,summary_target,fname,write_name)
# sys.exit()
#except:
# print("error")
#sys.exit()
#put into m_n list loop if adding to all concavity ksn calculations
m_n_list = [0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95]
for x in m_n_list:
#print(x)
x=str(x)
x = x.replace('.','_')
try:
os.rename(full_target+fname+'.bil',full_target+fname+x+'.bil')
os.rename(full_target+fname+'.hdr',full_target+fname+x+'.hdr')
except(OSError):
print("ERROR! Probably a missing tile. Location: %s%s"%(full_target,fname))
with open(full_target+fname+"_Chiculations.param",'wb') as file:
file.write('# This is a parameter file for the chi_mapping_tool \n')
file.write('# One day there will be documentation. \n')
file.write(" \n")
file.write('# These are parameters for the file i/o \n')
file.write("# IMPORTANT: You MUST make the write directory: the code will not work if it doens't exist. \n")
file.write('read path: %s \n'%(full_target))
file.write('write path: %s \n'%(full_target))
file.write('read fname: %s \n'%(fname+x))
file.write('write fname: %s \n'%(fname+x))
#file.write('read fname: %s \n'%(write_name))
#file.write('write fname: %s \n'%(write_name))
file.write("\n")
file.write('burn_raster_to_csv: True \n')
file.write('burn_raster_prefix: %s \n'%(fname+'_strain_ezz'))
#file.write('burn_raster_prefix: %s \n'%(fname+'_glaciated'))
#file.write('burn_data_csv_column_header: tertiary_burned_data \n')
file.write('burn_data_csv_column_header: second_inv \n')
file.write("\n")
file.write('secondary_burn_raster_to_csv: False \n')
file.write('secondary_burn_raster_prefix: %s \n'%(fname+'_distance_along'))
file.write('secondary_burn_data_csv_column_header: distance_along \n')
file.close()
#os.remove(full_target+fname+"_Chiculations.param")
#print("removed")
raster_burner_command = "/exports/csce/datastore/geos/users/s1134744/LSDTopoTools/Git_projects/LSDTT_Development/driver_functions_MuddChi2014/raster_burner.exe %s %s" %(full_target,fname+"_Chiculations.param")
#print(raster_burner_command)
try:
sub.call(raster_burner_command, shell = True)
#print(full_target+fname+x+'_MChiSegmented_burned.csv')
#target_dem = os.path.isfile(full_target+fname+x+'.bil')
#target_MChi = os.path.isfile(full_target+fname+x+'_MChiSegmented_burned.csv')
#target_burn = os.path.isfile(full_target+fname+'_to_burn.bil')
#print(target_dem,target_MChi,target_burn)
#if target_dem and target_MChi and target_burn:
# print("OK to run burn_raster")
except:
print("could not burn to csv!!!")
#undoing rename to allow renaming in next step
try:
os.rename(full_target+fname+x+'.bil',full_target+fname+'.bil')
os.rename(full_target+fname+x+'.hdr',full_target+fname+'.hdr')
except:
print("error fixing rename")
#clearing up
#os.remove(full_target+write_name+'_precipitation.bil')
#os.remove(full_target+write_name+'_precipitation.hdr')
#os.remove(full_target+write_name+'_geology.bil')
#os.remove(full_target+write_name+'_geology.hdr')