Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
_m = os.path.basename(trib.get_filepath('local_mustar')).split('.')
muin = _m[0] + in_id + '.' + _m[1]
muout = _m[0] + '_' + out_id + '.' + _m[1]
shutil.copyfile(os.path.join(trib.dir, muin),
os.path.join(main.dir, muout))
# sort flowlines descending
tfls.sort(key=lambda x: x.order, reverse=True)
# loop over tributaries and reproject to main glacier
for nr, tfl in enumerate(tfls):
# 1. Step: Change projection to the main glaciers grid
_line = salem.transform_geometry(tfl.line,
crs=trib.grid, to_crs=main.grid)
# 2. set new line
tfl.set_line(_line)
# 3. set map attributes
dx = [shpg.Point(tfl.line.coords[i]).distance(
shpg.Point(tfl.line.coords[i+1]))
for i, pt in enumerate(tfl.line.coords[:-1])] # get distance
# and check if equally spaced
if not np.allclose(dx, np.mean(dx), atol=1e-2):
raise RuntimeError('Flowline is not evenly spaced.')
tfl.dx = np.mean(dx).round(2)
tfl.map_dx = mfl.map_dx
tfl.dx_meter = tfl.map_dx * tfl.dx
def _extent_to_polygon(lon_ex, lat_ex, to_crs=None):
if lon_ex[0] == lon_ex[1] and lat_ex[0] == lat_ex[1]:
out = shpg.Point(lon_ex[0], lat_ex[0])
else:
x = [lon_ex[0], lon_ex[1], lon_ex[1], lon_ex[0], lon_ex[0]]
y = [lat_ex[0], lat_ex[0], lat_ex[1], lat_ex[1], lat_ex[0]]
out = shpg.Polygon(np.array((x, y)).T)
if to_crs is not None:
out = salem.transform_geometry(out, to_crs=to_crs)
return out
# open srtm tif-file:
dem = read_geotiff_dem(gdir)
if np.min(dem) == np.max(dem):
raise RuntimeError('({}) min equal max in the DEM.'
.format(gdir.rgi_id))
# Clip topography to 0 m a.s.l.
utils.clip_min(dem, 0, out=dem)
# Interpolate shape to a regular path
glacier_poly_hr = tolist(geometry)
for nr, poly in enumerate(glacier_poly_hr):
# transform geometry to map
_geometry = salem.transform_geometry(poly, to_crs=gdir.grid.proj)
glacier_poly_hr[nr] = _interp_polygon(_geometry, gdir.grid.dx)
glacier_poly_hr = shpg.MultiPolygon(glacier_poly_hr)
# Transform geometry into grid coordinates
# It has to be in pix center coordinates because of how skimage works
def proj(x, y):
grid = gdir.grid.center_grid
return grid.transform(x, y, crs=grid.proj)
glacier_poly_hr = shapely.ops.transform(proj, glacier_poly_hr)
# simple trick to correct invalid polys:
# http://stackoverflow.com/questions/20833344/
# fix-invalid-polygon-python-shapely
glacier_poly_hr = glacier_poly_hr.buffer(0)
# return list
tributaries = []
# loop over tributaries
for trib in candidates:
# skip self
if gdir.rgi_id == trib.rgi_id:
continue
# get tributary glacier downstream line and CRS
_dline = trib.read_pickle('downstream_line')['full_line']
_crs = trib.grid
# use salem to transform the grids
_trans_dline = salem.transform_geometry(_dline, crs=_crs, to_crs=crs)
# check for intersection, with a small buffer and add to list
if dline.intersects(_trans_dline.buffer(buffer)):
tributaries.append(trib)
return tributaries
mfls = main.read_pickle('model_flowlines')
# ------------------------------
# 0. create the new GlacierDirectory from main glaciers GeoDataFrame
# Should be passed along, if not download it
if glcdf is None:
glcdf = get_rgi_glacier_entities([main.rgi_id])
# Get index location of the specific glacier
idx = glcdf.loc[glcdf.RGIId == main.rgi_id].index
maindf = glcdf.loc[idx].copy()
# add tributary geometries to maindf
merged_geometry = maindf.loc[idx, 'geometry'].iloc[0]
for trib in tribs:
geom = trib.read_pickle('geometries')['polygon_hr']
geom = salem.transform_geometry(geom, crs=trib.grid)
merged_geometry = merged_geometry.union(geom)
# to get the center point, maximal extensions for DEM and single Polygon:
maindf.loc[idx, 'geometry'] = merged_geometry.convex_hull
# make some adjustments to the rgi dataframe
# 1. update central point of new glacier
# Could also use the mean of the Flowline centroids, to shift it more
# downstream. But this is more straight forward.
maindf.loc[idx, 'CenLon'] = maindf.loc[idx, 'geometry'].centroid.x
maindf.loc[idx, 'CenLat'] = maindf.loc[idx, 'geometry'].centroid.y
# 2. update names
maindf.loc[idx, 'RGIId'] += '_merged'
if maindf.loc[idx, 'Name'].iloc[0] is None:
maindf.loc[idx, 'Name'] = main.name + ' (merged)'
else: