@@ -168,62 +168,64 @@ def write_shape_file(fDict, shp_file, box=None, zero_first=False):
168
168
169
169
lats , lons = ut .get_lat_lon (ts_obj .metadata , geom_file = fDict ['Geometry' ], box = box )
170
170
171
- ###Loop over all datasets in context managers to skip close statements
172
- with h5py .File (fDict ['TimeSeries' ], 'r' ) as tsid :
173
- with h5py .File (fDict ['Coherence' ], 'r' ) as cohid :
174
- with h5py .File (fDict ['Velocity' ], 'r' ) as velid :
175
- with h5py .File (fDict ['Geometry' ], 'r' ) as geomid :
176
-
177
- length = box [3 ] - box [1 ]
178
- width = box [2 ] - box [0 ]
179
-
180
- #Start counter
181
- counter = 1
182
- prog_bar = ptime .progressBar (maxValue = nValid )
183
-
184
- #For each line
185
- for i in range (length ):
186
- line = i + box [1 ]
187
-
188
- # read data for the line
189
- ts = tsid ['timeseries' ][:, line , box [0 ]:box [2 ]].astype (np .float64 )
190
- if zero_first :
191
- ts -= np .tile (ts [0 , :], (ts .shape [0 ], 1 ))
192
-
193
- coh = cohid ['temporalCoherence' ][line , box [0 ]:box [2 ]].astype (np .float64 )
194
- vel = velid ['velocity' ][line , box [0 ]:box [2 ]].astype (np .float64 )
195
- vel_std = velid ['velocityStd' ][line , box [0 ]:box [2 ]].astype (np .float64 )
196
- hgt = geomid ['height' ][line , box [0 ]:box [2 ]].astype (np .float64 )
197
- lat = lats [i , :].astype (np .float64 )
198
- lon = lons [i , :].astype (np .float64 )
199
-
200
- for j in range (width ):
201
- if mask [i , j ] == 0 :
202
- continue
203
-
204
- #Create metadata dict
205
- rdict = { 'CODE' : hex (counter )[2 :].zfill (8 ),
206
- 'HEIGHT' : hgt [j ],
207
- 'H_STDEV' : 0. ,
208
- 'VEL' : vel [j ]* 1000 ,
209
- 'V_STDEV' : vel_std [j ]* 1000 ,
210
- 'COHERENCE' : coh [j ],
211
- 'EFF_AREA' : 1 }
212
-
213
- for ind , date in enumerate (ts_obj .dateList ):
214
- rdict [f'D{ date } ' ] = ts [ind , j ] * 1000
215
-
216
- #Create feature with definition
217
- feature = ogr .Feature (layerDefn )
218
- add_metadata (feature , [lon [j ], lat [j ]], rdict )
219
- layer .CreateFeature (feature )
220
- feature = None
221
-
222
- # update counter / progress bar
223
- counter += 1
224
- prog_bar .update (counter , every = 100 , suffix = f'line { counter } /{ nValid } ' )
225
- prog_bar .close ()
226
-
171
+ ### Use context managers to skip close statements
172
+ with (
173
+ h5py .File (fDict ["TimeSeries" ], "r" ) as tsid ,
174
+ h5py .File (fDict ["Coherence" ], "r" ) as cohid ,
175
+ h5py .File (fDict ["Velocity" ], "r" ) as velid ,
176
+ h5py .File (fDict ["Geometry" ], "r" ) as geomid ,
177
+ ):
178
+ length = box [3 ] - box [1 ]
179
+ width = box [2 ] - box [0 ]
180
+
181
+ # Start counter
182
+ counter = 1
183
+ prog_bar = ptime .progressBar (maxValue = nValid )
184
+
185
+ # For each line
186
+ for i in range (length ):
187
+ line = i + box [1 ]
188
+
189
+ # read data for the line
190
+ ts = tsid ["timeseries" ][:, line , box [0 ] : box [2 ]].astype (np .float64 )
191
+ if zero_first :
192
+ ts -= np .tile (ts [0 , :], (ts .shape [0 ], 1 ))
193
+
194
+ coh = cohid ["temporalCoherence" ][line , box [0 ] : box [2 ]].astype (np .float64 )
195
+ vel = velid ["velocity" ][line , box [0 ] : box [2 ]].astype (np .float64 )
196
+ vel_std = velid ["velocityStd" ][line , box [0 ] : box [2 ]].astype (np .float64 )
197
+ hgt = geomid ["height" ][line , box [0 ] : box [2 ]].astype (np .float64 )
198
+ lat = lats [i , :].astype (np .float64 )
199
+ lon = lons [i , :].astype (np .float64 )
200
+
201
+ for j in range (width ):
202
+ if mask [i , j ] == 0 :
203
+ continue
204
+
205
+ # Create metadata dict
206
+ rdict = {
207
+ "CODE" : hex (counter )[2 :].zfill (8 ),
208
+ "HEIGHT" : hgt [j ],
209
+ "H_STDEV" : 0.0 ,
210
+ "VEL" : vel [j ] * 1000 ,
211
+ "V_STDEV" : vel_std [j ] * 1000 ,
212
+ "COHERENCE" : coh [j ],
213
+ "EFF_AREA" : 1 ,
214
+ }
215
+
216
+ for ind , date in enumerate (ts_obj .dateList ):
217
+ rdict [f"D{ date } " ] = ts [ind , j ] * 1000
218
+
219
+ # Create feature with definition
220
+ feature = ogr .Feature (layerDefn )
221
+ add_metadata (feature , [lon [j ], lat [j ]], rdict )
222
+ layer .CreateFeature (feature )
223
+ feature = None
224
+
225
+ # update counter / progress bar
226
+ counter += 1
227
+ prog_bar .update (counter , every = 100 , suffix = f"line { counter } /{ nValid } " )
228
+ prog_bar .close ()
227
229
print (f'finished writing to file: { shp_file } ' )
228
230
return shp_file
229
231
0 commit comments