def transform(p1, p2, c1, c2, c3):
if PYPROJ:
if type(p1) is Proj and type(p2) is Proj:
if p1.srs != p2.srs:
return proj_transform(p1, p2, c1, c2, c3)
else:
return (c1, c2, c3)
elif type(p2) is TransverseMercator:
wgs84 = Proj(init="EPSG:4326")
if p1.srs != wgs84.srs:
t2, t1, t3 = proj_transform(p1, wgs84, c1, c2, c3)
else:
t1, t2, t3 = c2, c1, c3 # mind c2, c1 inversion
tm1, tm2 = p2.fromGeographic(t1, t2)
return (tm1, tm2, t3)
else:
if p1.spherical:
t1, t2 = p2.fromGeographic(c2, c1) # mind c2, c1 inversion
return (t1, t2, c3)
else:
return (c1, c2, c3)
python类transform()的实例源码
def convert_coordinates(coords, origin, wgs84, wrapped):
""" Convert coordinates from one crs to another """
if isinstance(coords, list) or isinstance(coords, tuple):
try:
if isinstance(coords[0], list) or isinstance(coords[0], tuple):
return [convert_coordinates(list(c), origin, wgs84, wrapped) for c in coords]
elif isinstance(coords[0], float):
c = list(transform(origin, wgs84, *coords))
if wrapped and c[0] < -170:
c[0] = c[0] + 360
return c
except IndexError:
pass
return None
def rasterizeShapes(pshapes, geodict, all_touched=True):
"""
Rastering a shape
Args:
pshapes: Sequence of orthographically projected shapes.
goedict: Mapio geodictionary.
all_touched: Turn pixel "on" if shape touches pixel, otherwise turn it
on if the center of the pixel is contained within the shape. Note
that the footprint of the shape is inflated and the amount of
inflation depends on the pixel size if all_touched=True.
Returns:
Rasterio grid.
"""
outshape = (geodict.ny, geodict.nx)
txmin = geodict.xmin - geodict.dx / 2.0
tymax = geodict.ymax + geodict.dy / 2.0
transform = Affine.from_gdal(
txmin, geodict.dx, 0.0, tymax, 0.0, -geodict.dy)
img = rasterio.features.rasterize(
pshapes, out_shape=outshape, fill=0.0, transform=transform,
all_touched=all_touched, default_value=1.0)
return img
def getClassBalance(pshapes, bounds, proj):
"""
Get native class balance of projected shapes, assuming a rectangular
bounding box.
Args:
pshapes: Sequence of projected shapely shapes.
bounds: Desired bounding box, in decimal degrees.
proj: PyProj object defining orthographic projection of shapes.
Returns:
Float fraction of hazard polygons (area of hazard polygons/total
area of bbox).
"""
xmin, ymin, xmax, ymax = bounds
bpoly = Polygon([(xmin, ymax),
(xmax, ymax),
(xmax, ymin),
(xmin, ymin)])
project = partial(
pyproj.transform,
pyproj.Proj(proj='latlong', datum='WGS84'),
proj)
bpolyproj = transform(project, bpoly)
totalarea = bpolyproj.area
polyarea = 0
for pshape in pshapes:
polyarea += pshape.area
return polyarea / totalarea
def get_angle_between_points(point1, point2, point3):
m_point1 = pyproj.transform(PROJ_WGS_84, PROJ_MERCATOR,
point1[LON], point1[LAT])
m_point2 = pyproj.transform(PROJ_WGS_84, PROJ_MERCATOR,
point2[LON], point2[LAT])
m_point3 = pyproj.transform(PROJ_WGS_84, PROJ_MERCATOR,
point3[LON], point3[LAT])
v1x = (m_point1[LON] - m_point2[LON]) # / COORDINATE_PRECISION
v1y = m_point1[LAT] - m_point2[LAT]
v2x = (m_point3[LON] - m_point2[LON]) # / COORDINATE_PRECISION
v2y = m_point3[LAT] - m_point2[LAT]
angle = (atan2(v2y, v2x) - atan2(v1y, v1x)) * 180. / pi
while angle < 0:
angle += 360.
return angle
def get_address_location(parcel_map_link):
"""
Parses the parcel map link and calculates coordinates from the extent.
An example link looks like this:
http://qpublic9.qpublic.net/qpmap4/map.php?county=la_orleans&parcel=41050873&extent=3667340+524208+3667804+524540&layers=parcels+aerials+roads+lakes
"""
o = urlparse(parcel_map_link)
query = parse_qs(o.query)
bbox = query['extent'][0].split(' ')
x1, y1, x2, y2 = [float(pt) for pt in bbox]
# get the midpoint of the extent
midpoint = [(x1 + x2) / 2, (y1 + y2) / 2]
# transform projected coordinates to latitude and longitude
in_proj = Proj(init='epsg:3452', preserve_units=True)
out_proj = Proj(init='epsg:4326')
return transform(in_proj, out_proj, midpoint[0], midpoint[1])
def transform(p1, p2, c1, c2, c3):
if PYPROJ:
if type(p1) is Proj and type(p2) is Proj:
if p1.srs != p2.srs:
return proj_transform(p1, p2, c1, c2, c3)
else:
return (c1, c2, c3)
elif type(p2) is TransverseMercator:
wgs84 = Proj(init="EPSG:4326")
if p1.srs != wgs84.srs:
t2, t1, t3 = proj_transform(p1, wgs84, c1, c2, c3)
else:
t1, t2, t3 = c2, c1, c3 # mind c2, c1 inversion
tm1, tm2 = p2.fromGeographic(t1, t2)
return (tm1, tm2, t3)
else:
if p1.spherical:
t1, t2 = p2.fromGeographic(c2, c1) # mind c2, c1 inversion
return (t1, t2, c3)
else:
return (c1, c2, c3)
def latlon_to_ij(self, lat, lon):
"""
Convert latitude and longitude into grid coordinates.
If this is a child domain, it asks it's parent to do the projectiona and then
remaps it into its own coordinate system via parent_start and cell size ratio.
:param lat: latitude
:param lon: longitude
:return: the i, j position in grid coordinates
"""
if self.top_level:
proj_j, proj_i = pyproj.transform(self.latlon_sphere, self.lambert_grid,
lon, lat)
return ((proj_i - self.offset_i) / self.cell_size[0],
(proj_j - self.offset_j) / self.cell_size[1])
else:
pi, pj = self.parent.latlon_to_ij(lat, lon)
pcsr, ps = self.parent_cell_size_ratio, self.parent_start
delta = (pcsr - 1) / 2
return ((pi - ps[0] + 1.) * pcsr + delta,
(pj - ps[1] + 1.) * pcsr + delta)
def txncsptolatlon (northing, easting) :
'''
Sweetwater
Convert texas state plane coordinates in feet to
geographic coordinates, WGS84.
'''
# Texas NC state plane feet Zone 4202
sp = Proj (init='epsg:32038')
# WGS84, geographic
wgs = Proj (init='epsg:4326', proj='latlong')
# Texas SP coordinates: survey foot is 1200/3937 meters
lon, lat = transform (sp, wgs, easting * 0.30480060960121924, northing * 0.30480060960121924)
return lat, lon
def utmcsptolatlon (northing, easting) :
'''
Mount Saint Helens
Convert UTM to
geographic coordinates, WGS84.
'''
# UTM
utmc = Proj (proj='utm', zone=UTM, ellps='WGS84')
# WGS84, geographic
wgs = Proj (init='epsg:4326', proj='latlong')
#
lon, lat = transform (utmc, wgs, easting, northing)
return lat, lon
def __getitem__(self, geometry):
if isinstance(geometry, BaseGeometry) or getattr(geometry, "__geo_interface__", None) is not None:
image = GeoImage.__getitem__(self, geometry)
return image
else:
result = DaskImage.__getitem__(self, geometry)
image = super(GeoDaskWrapper, self.__class__).__new__(self.__class__,
result.dask, result.name, result.chunks,
result.dtype, result.shape)
if all([isinstance(e, slice) for e in geometry]) and len(geometry) == len(self.shape):
xmin, ymin, xmax, ymax = geometry[2].start, geometry[1].start, geometry[2].stop, geometry[1].stop
xmin = 0 if xmin is None else xmin
ymin = 0 if ymin is None else ymin
xmax = self.shape[2] if xmax is None else xmax
ymax = self.shape[1] if ymax is None else ymax
g = ops.transform(self.__geo_transform__.fwd, box(xmin, ymin, xmax, ymax))
image.__geo_interface__ = mapping(g)
image.__geo_transform__ = self.__geo_transform__ + (xmin, ymin)
else:
image.__geo_interface__ = self.__geo_interface__
image.__geo_transform__ = self.__geo_transform__
return image
def transform_coordinates(self, to_projection):
try:
from pyproj import Proj, transform
inProj = Proj(init=self.metadata['projection'])
outProj = Proj(init=to_projection)
x, y = transform(inProj, outProj, self.metadata['x'],
self.metadata['y'])
self.metadata['x'] = x
self.metadata['y'] = y
self.metadata['projection'] = to_projection
except ImportError:
raise ImportError(
'The module pyproj could not be imported. Please '
'install through:'
'>>> pip install pyproj'
'or ... conda install pyproj')
def xyz2geodetic(x, y, z):
"""
Convert ECEF *x*, *y*, and *z* (all in [m]) to geodetic
coordinates (using the WGS84 ellipsoid). Return lat [deg], lon
[deg], and alt [m]. Multiple coordinates are acceptable, i.e.,
*x*, *y*, and *z* may be equal length containers.
"""
lon, lat, alt = pyproj.transform(ECEF, LLA, x, y, z, radians=False)
if isinstance(lon, Iterable):
lon = [x - 360 if x > 180 else x for x in lon]
else:
if lon > 180:
lon -= 360
return lat, lon, alt
def pixel_to_lat_lon_web_mercator(raster_dataset, col, row):
"""Convert a pixel on the raster_dataset to web mercator (epsg:3857)."""
spatial_reference = osr.SpatialReference()
spatial_reference.ImportFromWkt(raster_dataset.GetProjection())
ds_spatial_reference_proj_string = spatial_reference.ExportToProj4()
in_proj = Proj(ds_spatial_reference_proj_string)
out_proj = Proj(init='epsg:3857')
geo_transform = raster_dataset.GetGeoTransform()
ulon = col * geo_transform[1] + geo_transform[0]
ulat = row * geo_transform[5] + geo_transform[3]
x2, y2 = transform(in_proj, out_proj, ulon, ulat)
x2, y2 = out_proj(x2, y2, inverse=True)
return x2, y2
def test_wrap_coordinates(coords, origin, wgs84):
""" Test whether coordinates wrap around the antimeridian in wgs84 """
lon_under_minus_170 = False
lon_over_plus_170 = False
if isinstance(coords[0], list):
for c in coords[0]:
c = list(transform(origin, wgs84, *c))
if c[0] < -170:
lon_under_minus_170 = True
elif c[0] > 170:
lon_over_plus_170 = True
else:
return False
return lon_under_minus_170 and lon_over_plus_170
def latlon2pixel(lat, lon, input_raster='', targetsr='', geom_transform=''):
# type: (object, object, object, object, object) -> object
sourcesr = osr.SpatialReference()
sourcesr.ImportFromEPSG(4326)
geom = ogr.Geometry(ogr.wkbPoint)
geom.AddPoint(lon, lat)
if targetsr == '':
src_raster = gdal.Open(input_raster)
targetsr = osr.SpatialReference()
targetsr.ImportFromWkt(src_raster.GetProjectionRef())
coord_trans = osr.CoordinateTransformation(sourcesr, targetsr)
if geom_transform == '':
src_raster = gdal.Open(input_raster)
transform = src_raster.GetGeoTransform()
else:
transform = geom_transform
x_origin = transform[0]
# print(x_origin)
y_origin = transform[3]
# print(y_origin)
pixel_width = transform[1]
# print(pixel_width)
pixel_height = transform[5]
# print(pixel_height)
geom.Transform(coord_trans)
# print(geom.GetPoint())
x_pix = (geom.GetPoint()[0] - x_origin) / pixel_width
y_pix = (geom.GetPoint()[1] - y_origin) / pixel_height
return (x_pix, y_pix)
def osgb_point(easting, northing):
return "[%.5f,%.5f]" % pyproj.transform(osgb36, wgs84, easting, northing)
def osgb_point(easting, northing):
return "[%.5f,%.5f]" % pyproj.transform(osgb36, wgs84, easting, northing)
# load map of SNAC LA codes to the local-authority register
def getProjectedShapes(shapes, xmin, xmax, ymin, ymax):
"""
Take a sequence of geographic shapes and project them to a bounds-centered
orthographic projection.
Args:
shapes: Sequence of shapes, as read in by fiona.collection().
xmin: Eastern boundary of all shapes.
xmax: Western boundary of all shapes.
ymin: Southern boundary of all shapes.
ymax: Northern boundary of all shapes.
Returns:
Tuple of
- Input sequence of shapes, projected to orthographic
- PyProj projection object used to transform input shapes
"""
latmiddle = ymin + (ymax - ymin) / 2.0
lonmiddle = xmin + (xmax - xmin) / 2.0
projstr = ('+proj=ortho +datum=WGS84 +lat_0=%.4f +lon_0=%.4f '
'+x_0=0.0 +y_0=0.0' % (latmiddle, lonmiddle))
proj = pyproj.Proj(projparams=projstr)
project = partial(
pyproj.transform,
pyproj.Proj(proj='latlong', datum='WGS84'),
proj)
pshapes = []
for tshape in shapes:
if tshape['geometry']['type'] == 'Polygon':
pshapegeo = shape(tshape['geometry'])
else:
pshapegeo = shape(tshape['geometry'])
pshape = transform(project, pshapegeo)
pshapes.append(pshape) # assuming here that these are simple polygons
return (pshapes, proj)
def projectBack(points, proj):
"""
Project a 2D array of XY points from orthographic projection to decimal
degrees.
Args:
points: 2D numpy array of XY points in orthographic projection.
proj: PyProj object defining projection.
Returns:
2D numpy array of Lon/Lat coordinates.
"""
mpoints = MultiPoint(points)
project = partial(
pyproj.transform,
proj,
pyproj.Proj(proj='latlong', datum='WGS84'))
gmpoints = transform(project, mpoints)
coords = []
for point in gmpoints.geoms:
x, y = point.coords[0]
coords.append((x, y))
coords = np.array(coords)
return coords
def projectBack(points, proj):
"""
Project a 2D array of XY points from orthographic projection to decimal
degrees.
Args:
points: 2D numpy array of XY points in orthographic projection.
proj: PyProj object defining projection.
Returns:
2D numpy array of Lon/Lat coordinates.
"""
mpoints = MultiPoint(points)
project = partial(
pyproj.transform,
proj,
pyproj.Proj(proj='latlong', datum='WGS84'))
gmpoints = transform(project, mpoints)
coords = []
for point in gmpoints.geoms:
x, y = point.coords[0]
coords.append((x, y))
coords = np.array(coords)
return coords
def area_from_lon_lat_poly(geometry):
"""
Compute the area in km^2 of a shapely geometry, whose points are in longitude and latitude.
Parameters
----------
geometry: shapely geometry
Points must be in longitude and latitude.
Returns
-------
area: float
Area in km^2.
"""
import pyproj
from shapely.ops import transform
from functools import partial
project = partial(
pyproj.transform,
pyproj.Proj(init='epsg:4326'), # Source: Lon-Lat
pyproj.Proj(proj='aea')) # Target: Albers Equal Area Conical https://en.wikipedia.org/wiki/Albers_projection
new_geometry = transform(project, geometry)
#default area is in m^2
return new_geometry.area/1e6
def get_proj(epsg):
"""Returns a pyproj partial representing a projedction from WGS84 to
a given projection.
Args:
epsg: EPSG code for the target projection.
"""
project = partial(
pyproj.transform,
pyproj.Proj(init='EPSG:4326'),
pyproj.Proj(init='EPSG:{epsg}'.format(epsg=epsg)))
def ny_state_coords_to_lat_lng(xcoord, ycoord):
return pyproj.transform(ny_state_plane, wgs84, xcoord, ycoord)
def ij_to_latlon(self, i, j):
if self.top_level:
lon, lat = pyproj.transform(self.lambert_grid, self.latlon_sphere,
j * self.cell_size[1] + self.offset_j,
i * self.cell_size[0] + self.offset_i)
return lat, lon
else:
pcsr, ps = self.parent_cell_size_ratio, self.parent_start
delta = (pcsr - 1) / 2
return self.parent.ij_to_latlon((i-delta)/pcsr+ps[0]-1., (j-delta)/pcsr+ps[1]-1.)
def convert_to_latlon(geoimg, lines):
srs = osr.SpatialReference(geoimg.srs()).ExportToProj4()
projin = Proj(srs)
projout = Proj(init='epsg:4326')
newlines = []
for line in lines:
l = []
for point in line:
pt = transform(projin, projout, point[0], point[1])
l.append(pt)
newlines.append(l)
return antimeridian_linesplit(newlines)
def bbox_3857(self):
inProj = Proj(init='epsg:4326')
outProj = Proj(init='epsg:3857')
sw = transform(inProj, outProj, self.bbox_x0, self.bbox_y0)
ne = transform(inProj, outProj, self.bbox_x1, self.bbox_y1)
return [sw[0], sw[1], ne[0], ne[1]]
def project(geom, from_proj=None, to_proj=None):
"""
Project a ``shapely`` geometry, and returns a new geometry of the same type from the transformed coordinates.
Default input projection is `WGS84 <https://en.wikipedia.org/wiki/World_Geodetic_System>`_, default output projection is `Mollweide <http://spatialreference.org/ref/esri/54009/>`_.
Inputs:
*geom*: A ``shapely`` geometry.
*from_proj*: A ``PROJ4`` string. Optional.
*to_proj*: A ``PROJ4`` string. Optional.
Returns:
A ``shapely`` geometry.
"""
from_proj = wgs84(from_proj)
if to_proj is None:
to_proj = MOLLWEIDE
else:
to_proj = wgs84(to_proj)
to_proj, from_proj = pyproj.Proj(to_proj), pyproj.Proj(from_proj)
if ((to_proj == from_proj) or
(to_proj.is_latlong() and from_proj.is_latlong())):
return geom
projection_func = partial(pyproj.transform, from_proj, to_proj)
return transform(projection_func, geom)
def convert_to_ecef(x, y, z, epsg_input):
inp = pyproj.Proj(init='epsg:{0}'.format(epsg_input))
outp = pyproj.Proj(init='epsg:4978') # ECEF
return pyproj.transform(inp, outp, x, y, z)
def _tile_coords(self, bounds):
""" Convert tile coords mins/maxs to lng/lat bounds """
tfm = partial(pyproj.transform,
pyproj.Proj(init="epsg:3857"),
pyproj.Proj(init="epsg:4326"))
bounds = ops.transform(tfm, box(*bounds)).bounds
params = list(bounds) + [[self.zoom_level]]
tile_coords = [(tile.x, tile.y) for tile in mercantile.tiles(*params)]
xtiles, ytiles = zip(*tile_coords)
minx = min(xtiles)
maxx = max(xtiles)
miny = min(ytiles)
maxy = max(ytiles)
return minx, miny, maxx, maxy
