Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def test_transform_direction__string():
forward_transformer = Transformer.from_crs(4326, 3857)
inverse_transformer = Transformer.from_crs(3857, 4326)
assert inverse_transformer.transform(
-33, 24, direction="INVERSE"
) == forward_transformer.transform(-33, 24, direction="FORWARD")
ident_transformer = Transformer.from_crs(4326, 3857)
ident_transformer.transform(-33, 24, direction="IDENT") == (-33, 24)
def project(shape, source, target):
"""Projects a geometry from one coordinate system into another.
Args:
shape: the geometry to project.
source: the source EPSG spatial reference system identifier.
target: the target EPSG spatial reference system identifier.
Returns:
The projected geometry in the target coordinate system.
"""
transformer = pyproj.Transformer.from_crs(source, target)
return shapely.ops.transform(transformer.transform, shape)
def _tile_coords(self, bounds):
""" convert mercator bbox to tile index limits """
tfm = pyproj.Transformer.from_crs(3857, 4326, always_xy=True)
bounds = ops.transform(tfm.transform, box(*bounds)).bounds
# because tiles have a common corner, the tiles that cover a
# given tile includes the adjacent neighbors.
# https://github.com/mapbox/mercantile/issues/84#issuecomment-413113791
west, south, east, north = bounds
epsilon = 1.0e-10
if east != west and north != south:
# 2D bbox
# shrink the bounds a small amount so that
# shapes/tiles round trip.
west += epsilon
south += epsilon
east -= epsilon
north -= epsilon
if isinstance(other, Airspace):
other = other.flatten()
if isinstance(other, Polygon):
bounds = other.bounds
projection = pyproj.Proj(
proj="aea", # equivalent projection
lat_1=bounds[1],
lat_2=bounds[3],
lat_0=(bounds[1] + bounds[3]) / 2,
lon_0=(bounds[0] + bounds[2]) / 2,
)
transformer = pyproj.Transformer.from_proj(
pyproj.Proj("epsg:4326"), projection, always_xy=True
)
projected_shape = transform(transformer.transform, other)
self_xy = self.compute_xy(projection)
return self.assign(
**{
column_name: list(
projected_shape.exterior.distance(p)
* (-1 if projected_shape.contains(p) else 1)
for p in MultiPoint(
list(zip(self_xy.data.x, self_xy.data.y))
)
)
}
"""Snowfall analysis maps."""
import datetime
import numpy as np
import pandas as pd
from pyproj import Transformer
from geopandas import read_postgis, GeoDataFrame
from shapely.geometry import Polygon, Point
from scipy.interpolate import Rbf
from scipy.ndimage import zoom
from pyiem import reference
from pyiem.plot import MapPlot, nwssnow
from pyiem.util import get_autoplot_context, get_dbconn
T4326_2163 = Transformer.from_proj(4326, 2163, always_xy=True)
T2163_4326 = Transformer.from_proj(2163, 4326, always_xy=True)
USEME = "used_for_analysis"
SCIPY = (
"https://docs.scipy.org/doc/scipy/reference/generated/"
"scipy.interpolate.Rbf.html"
)
PDICT = {"cwa": "Plot by NWS Forecast Office", "state": "Plot by State"}
PDICT2 = {
"multiquadric": "multiquadraic sqrt((r/self.epsilon)**2 + 1)",
"inverse": "inverse 1.0/sqrt((r/self.epsilon)**2 + 1)",
"gaussian": "gaussian exp(-(r/self.epsilon)**2)",
"linear": "linear r",
"cubic": "cubic r**3",
"quintic": "quintic r**5",
"thin_plate": "thin_plate r**2 * log(r)",
}
PDICT3 = {
def draw_ruler(raster, left, top, right, bot, step, crs, font, color, w, dh):
img = Image.fromarray(raster)
d = ImageDraw.Draw(img)
def f(x, dir=0):
v = raster.shape[1-dir]
if abs(x)>1: x = int(x)
else: x = int(v * x)
return x if x>=0 else v+x
left, top, right, bot = f(left), f(top,1), f(right), f(bot,1)
d.rectangle([left, top, right, bot], outline=color, width=w)
pts = np.array([(left,top),(right,top),(right,bot),(left,bot)])
prj1, prj2 = pyproj.CRS(raster.crs), pyproj.CRS(crs)
ct = pyproj.Transformer.from_crs(prj1, prj2, always_xy=True)
pts = np.dot(raster.mat[:,1:], pts.T) + raster.mat[:,:1]
xs, ys = ct.transform(*pts)
a, b = int(xs[0]//step), int(xs[1]//step)
tab = [i*step for i in range(a+1, b+1)]
nps = ct.transform(tab, np.linspace(ys[0], ys[1], len(tab)), direction='INVERSE')
nps = np.dot(np.linalg.inv(raster.mat[:,1:]), nps - raster.mat[:,:1])
font = ImageFont.truetype(*font)
for x,e in zip(nps[0],tab):
d.line([int(x), top, int(x), top-dh], color, w)
tw, th = d.textsize('%dยฐE'%e, font)
d.text((int(x-tw//2), top-dh-th*1.2), '%dยฐE'%e, font=font, fill=color)
a, b = int(xs[3]//step), int(xs[2]//step)
tab = [i*step for i in range(a+1, b+1)]
nps = ct.transform(tab, np.linspace(ys[3], ys[2], len(tab)), direction='INVERSE')
def execute(cls, args, env):
tc_ids = DBSession.query(ResourceTileCache.resource_id).filter(
ResourceTileCache.enabled,
ResourceTileCache.seed_z != None # NOQA: E711
).all()
# TODO: Add arbitrary SRS support
srs_tr = Transformer.from_crs(4326, 3857, always_xy=True)
for tc_id in tc_ids:
tc = ResourceTileCache.filter_by(resource_id=tc_id).one()
rend_res = tc.resource
data_res = rend_res.parent
srs = data_res.srs
# TODO: Add arbitrary SRS support
extent_4326 = data_res.extent
extent = srs_tr.transform(extent_4326['minLon'], extent_4326['minLat']) + \
srs_tr.transform(extent_4326['maxLon'], extent_4326['maxLat'])
rlevel = list()
rcount = 0
if tolerance < 0:
raise ValueError("tolerance must be a positive float")
if df is not None and isinstance(lat, str) and isinstance(lon, str):
lat, lon = df[lat], df[lon]
if df is not None and lat is not None and lon is not None:
projection = pyproj.Proj(
proj="lcc",
ellps="WGS84",
lat_1=lat.min(),
lat_2=lat.max(),
lat_0=lat.mean(),
lon_0=lon.mean(),
)
transformer = pyproj.Transformer.from_proj(
pyproj.Proj("epsg:4326"), projection, always_xy=True
)
x, y = transformer.transform(lon.values, lat.values,)
else:
if df is not None:
x, y = df[x].values, df[y].values
x, y = np.array(x), np.array(y)
if z is not None:
if df is not None:
z = df[z].values
z = z_factor * np.array(z)
mask = np.ones(len(x), dtype=bool)
if z is None:
_douglas_peucker_rec(x, y, mask, tolerance)
return None
if isinstance(projection, crs.Projection):
projection = pyproj.Proj(projection.proj4_init)
if projection is None:
bounds = self.bounds
projection = pyproj.Proj(
proj="aea", # equivalent projection
lat_1=bounds[1],
lat_2=bounds[3],
lat_0=(bounds[1] + bounds[3]) / 2,
lon_0=(bounds[0] + bounds[2]) / 2,
)
transformer = pyproj.Transformer.from_proj(
pyproj.Proj("epsg:4326"), projection, always_xy=True
)
projected_shape = transform(transformer.transform, self.shape,)
if not projected_shape.is_valid:
warnings.warn("The chosen projection is invalid for current shape")
return projected_shape
"""Snowfall analysis maps."""
import datetime
import numpy as np
import pandas as pd
from pyproj import Transformer
from geopandas import read_postgis, GeoDataFrame
from shapely.geometry import Polygon, Point
from scipy.interpolate import Rbf
from scipy.ndimage import zoom
from pyiem import reference
from pyiem.plot import MapPlot, nwssnow
from pyiem.util import get_autoplot_context, get_dbconn
T4326_2163 = Transformer.from_proj(4326, 2163, always_xy=True)
T2163_4326 = Transformer.from_proj(2163, 4326, always_xy=True)
USEME = "used_for_analysis"
SCIPY = (
"https://docs.scipy.org/doc/scipy/reference/generated/"
"scipy.interpolate.Rbf.html"
)
PDICT = {"cwa": "Plot by NWS Forecast Office", "state": "Plot by State"}
PDICT2 = {
"multiquadric": "multiquadraic sqrt((r/self.epsilon)**2 + 1)",
"inverse": "inverse 1.0/sqrt((r/self.epsilon)**2 + 1)",
"gaussian": "gaussian exp(-(r/self.epsilon)**2)",
"linear": "linear r",
"cubic": "cubic r**3",
"quintic": "quintic r**5",
"thin_plate": "thin_plate r**2 * log(r)",
}