Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def transform(self, x, y, z=None, crs=wgs84, nearest=False, maskout=False):
"""Converts any coordinates into the local grid.
Parameters
----------
x : ndarray
the grid coordinates of the point(s) you want to convert
y : ndarray
the grid coordinates of the point(s) you want to convert
z : None
ignored (but necessary since some shapes have a z dimension)
crs : crs
reference system of x, y. Could be a pyproj.Proj instance or a
Grid instance. In the latter case (x, y) are actually (i, j).
(Default: lonlat in wgs84).
nearest : bool
set to True if you wish to return the closest i, j coordinates
-------
the topography if needed (bonus)
"""
if topo is None:
self._shading_base()
return
kwargs.setdefault('interp', 'spline')
if isinstance(topo, six.string_types):
_, ext = os.path.splitext(topo)
if ext.lower() == '.tif':
g = GeoTiff(topo)
# Spare memory
ex = self.grid.extent_in_crs(crs=wgs84) # l, r, b, t
g.set_subset(corners=((ex[0], ex[2]), (ex[1], ex[3])),
crs=wgs84, margin=10)
z = g.get_vardata()
z[z < -999] = 0
z = self.grid.map_gridded_data(z, g.grid, **kwargs)
else:
raise ValueError('File extension not recognised: {}'
.format(ext))
else:
z = self._check_data(topo, crs=crs, **kwargs)
# Gradient in m m-1
ddx = self.grid.dx
ddy = self.grid.dy
if gis.proj_is_latlong(self.grid.proj):
# we make a coarse approx of the avg dx on a sphere
if gdir.is_tidewater and cfg.PARAMS['clip_tidewater_border']:
border = 10
# Corners, incl. a buffer of N pix
ulx = np.min(xx) - border * dx
lrx = np.max(xx) + border * dx
uly = np.max(yy) + border * dx
lry = np.min(yy) - border * dx
# n pixels
nx = np.int((lrx - ulx) / dx)
ny = np.int((uly - lry) / dx)
# Back to lon, lat for DEM download/preparation
tmp_grid = salem.Grid(proj=proj_out, nxny=(nx, ny), x0y0=(ulx, uly),
dxdy=(dx, -dx), pixel_ref='corner')
minlon, maxlon, minlat, maxlat = tmp_grid.extent_in_crs(crs=salem.wgs84)
# Open DEM
source = entity.DEM_SOURCE if hasattr(entity, 'DEM_SOURCE') else None
if not is_dem_source_available(source,
(minlon, maxlon),
(minlat, maxlat)):
raise InvalidDEMError('Source: {} not available for glacier {}'
.format(source, gdir.rgi_id))
dem_list, dem_source = get_topo_file((minlon, maxlon), (minlat, maxlat),
rgi_region=gdir.rgi_region,
rgi_subregion=gdir.rgi_subregion,
dx_meter=dx,
source=source)
log.debug('(%s) DEM source: %s', gdir.rgi_id, dem_source)
log.debug('(%s) N DEM Files: %s', gdir.rgi_id, len(dem_list))
# OK, get it
lon = ds.variables[x][:]
lat = ds.variables[y][:]
# double check for dubious variables
if not utils.str_in_list([x], utils.valid_names['lon_var']) or \
not utils.str_in_list([y], utils.valid_names['lat_var']):
# name not usual. see if at least the range follows some conv
if (np.max(np.abs(lon)) > 360.1) or (np.max(np.abs(lat)) > 90.1):
return None
# Make the grid
dx = lon[1]-lon[0]
dy = lat[1]-lat[0]
args = dict(nxny=(lon.shape[0], lat.shape[0]), proj=wgs84, dxdy=(dx, dy),
x0y0=(lon[0], lat[0]))
return gis.Grid(**args)
def is_rotated_proj_working():
import pyproj
srs = ('+ellps=WGS84 +proj=ob_tran +o_proj=latlon '
'+to_meter=0.0174532925199433 +o_lon_p=0.0 +o_lat_p=80.5 '
'+lon_0=357.5 +no_defs')
p1 = pyproj.Proj(srs)
p2 = wgs84
return np.isclose(transform_proj(p1, p2, -20, -9),
[-22.243473889042903, -0.06328365194179102],
atol=1e-5).all()
def googlestatic_mercator_grid(center_ll=None, nx=640, ny=640, zoom=12):
"""Mercator map centered on a specified point as seen by google API"""
# Make a local proj
lon, lat = center_ll
proj_params = dict(proj='merc', datum='WGS84')
projloc = pyproj.Proj(proj_params)
# Meter per pixel
mpix = (2 * np.pi * google_earth_radius) / google_pix / (2**zoom)
xx = nx * mpix
yy = ny * mpix
e, n = pyproj.transform(wgs84, projloc, lon, lat)
corner = (-xx / 2. + e, yy / 2. + n)
dxdy = (xx / nx, - yy / ny)
return Grid(proj=projloc, corner=corner, nxny=(nx, ny), dxdy=dxdy,
pixel_ref='corner')