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类Proj()的实例源码
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 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 __init__(self):
self.path = os.path.join(os.path.dirname(__file__),'data')
self.stat_data = {}
self.header_conversion = {}
self.stat_data, self.header_conversion = \
json.load(open(os.path.join(self.path,'stat-data.json')))
geometries = json.load(open(os.path.join(self.path,'geometries.json')))
crs = json.load(open(os.path.join(self.path,'crs.json')))
self.proj = pyproj.Proj(crs)
self.recs = [{
'shp':shapely.geometry.asShape(geometry),
'area': area_id,
} for area_id, geometry in geometries.items()]
for r in self.recs:
bounds = r['shp'].bounds
r['key'] = bounds[2] + bounds[3]
r['bounds'] = bounds
self.recs.sort(key=lambda r:r['key'])
self.rec_keys = [r['key'] for r in self.recs]
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 geod2utm (zn, datum, lat, lon, elev) :
''' Convert geodetic coordinates to UTM '''
if zn == None :
zn = lon2zone (lon)
p = Proj (proj='utm', zone=zn, ellps=datum)
X, Y = p (lon, lat)
# Return Y, X, Z
return Y, X, elev
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 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 init_proj(self, latlon):
self.proj = Proj(proj='utm',zone=calculate_utm_zone(*latlon)[0],ellps='WGS84')
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 condition_df(self):
"""
Do any initial data conditioning that may be required.
"""
logging.info('Ensure that columns that are supposed to be numeric are numeric')
self.df[SET_GHI] = pd.to_numeric(self.df[SET_GHI], errors='coerce')
self.df[SET_WINDVEL] = pd.to_numeric(self.df[SET_WINDVEL], errors='coerce')
self.df[SET_NIGHT_LIGHTS] = pd.to_numeric(self.df[SET_NIGHT_LIGHTS], errors='coerce')
self.df[SET_ELEVATION] = pd.to_numeric(self.df[SET_ELEVATION], errors='coerce')
self.df[SET_SLOPE] = pd.to_numeric(self.df[SET_SLOPE], errors='coerce')
self.df[SET_LAND_COVER] = pd.to_numeric(self.df[SET_LAND_COVER], errors='coerce')
self.df[SET_GRID_DIST_CURRENT] = pd.to_numeric(self.df[SET_GRID_DIST_CURRENT], errors='coerce')
self.df[SET_GRID_DIST_PLANNED] = pd.to_numeric(self.df[SET_GRID_DIST_PLANNED], errors='coerce')
self.df[SET_SUBSTATION_DIST] = pd.to_numeric(self.df[SET_SUBSTATION_DIST], errors='coerce')
self.df[SET_ROAD_DIST] = pd.to_numeric(self.df[SET_ROAD_DIST], errors='coerce')
self.df[SET_HYDRO_DIST] = pd.to_numeric(self.df[SET_HYDRO_DIST], errors='coerce')
self.df[SET_HYDRO] = pd.to_numeric(self.df[SET_HYDRO], errors='coerce')
self.df[SET_SOLAR_RESTRICTION] = pd.to_numeric(self.df[SET_SOLAR_RESTRICTION], errors='coerce')
logging.info('Replace null values with zero')
self.df.fillna(0, inplace=True)
logging.info('Sort by country, Y and X')
self.df.sort_values(by=[SET_COUNTRY, SET_Y, SET_X], inplace=True)
logging.info('Add columns with location in degrees')
project = Proj('+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')
def get_x(row):
x, y = project(row[SET_X] * 1000, row[SET_Y] * 1000, inverse=True)
return x
def get_y(row):
x, y = project(row[SET_X] * 1000, row[SET_Y] * 1000, inverse=True)
return y
self.df[SET_X_DEG] = self.df.apply(get_x, axis=1)
self.df[SET_Y_DEG] = self.df.apply(get_y, axis=1)
def condition_df(self):
"""
Do any initial data conditioning that may be required.
"""
logging.info('Ensure that columns that are supposed to be numeric are numeric')
self.df[SET_GHI] = pd.to_numeric(self.df[SET_GHI], errors='coerce')
self.df[SET_WINDVEL] = pd.to_numeric(self.df[SET_WINDVEL], errors='coerce')
self.df[SET_NIGHT_LIGHTS] = pd.to_numeric(self.df[SET_NIGHT_LIGHTS], errors='coerce')
self.df[SET_ELEVATION] = pd.to_numeric(self.df[SET_ELEVATION], errors='coerce')
self.df[SET_SLOPE] = pd.to_numeric(self.df[SET_SLOPE], errors='coerce')
self.df[SET_LAND_COVER] = pd.to_numeric(self.df[SET_LAND_COVER], errors='coerce')
self.df[SET_GRID_DIST_CURRENT] = pd.to_numeric(self.df[SET_GRID_DIST_CURRENT], errors='coerce')
self.df[SET_GRID_DIST_PLANNED] = pd.to_numeric(self.df[SET_GRID_DIST_PLANNED], errors='coerce')
self.df[SET_SUBSTATION_DIST] = pd.to_numeric(self.df[SET_SUBSTATION_DIST], errors='coerce')
self.df[SET_ROAD_DIST] = pd.to_numeric(self.df[SET_ROAD_DIST], errors='coerce')
self.df[SET_HYDRO_DIST] = pd.to_numeric(self.df[SET_HYDRO_DIST], errors='coerce')
self.df[SET_HYDRO] = pd.to_numeric(self.df[SET_HYDRO], errors='coerce')
self.df[SET_SOLAR_RESTRICTION] = pd.to_numeric(self.df[SET_SOLAR_RESTRICTION], errors='coerce')
logging.info('Replace null values with zero')
self.df.fillna(0, inplace=True)
logging.info('Sort by country, Y and X')
self.df.sort_values(by=[SET_COUNTRY, SET_Y, SET_X], inplace=True)
logging.info('Add columns with location in degrees')
project = Proj('+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')
def get_x(row):
x, y = project(row[SET_X] * 1000, row[SET_Y] * 1000, inverse=True)
return x
def get_y(row):
x, y = project(row[SET_X] * 1000, row[SET_Y] * 1000, inverse=True)
return y
self.df[SET_X_DEG] = self.df.apply(get_x, axis=1)
self.df[SET_Y_DEG] = self.df.apply(get_y, axis=1)
def georeference(self, scene, center):
if "latitude" not in scene and "longitude" not in scene:
if type(self.pScene) is TransverseMercator:
scene['latitude'] = self.pScene.lat
scene['longitude'] = self.pScene.lon
scene['altitude'] = 0
elif type(self.pScene) is not None:
wgs84 = Proj(init="EPSG:4326")
latlon = transform(self.pScene, wgs84, center[0], center[1], center[2])
scene['longitude'] = latlon[0]
scene['latitude'] = latlon[1]
scene['altitude'] = latlon[2]
def entities(self, name, scene=None):
"""
Iterates over all DXF entities according to the options set by user.
"""
if scene is None:
scene = bpy.context.scene
self.current_scene = scene
if self.recenter:
self.objects_before += scene.objects[:]
if self.combination == BY_BLOCK:
self.combined_objects((en for en in self.dwg.modelspace()), scene)
elif self.combination != SEPARATED:
self.combined_objects((en for en in self.dwg.modelspace() if is_.combined_entity(en)), scene)
self.separated_entities((en for en in self.dwg.modelspace() if is_.separated_entity(en)), scene)
else:
self.separated_entities((en for en in self.dwg.modelspace() if en.dxftype != "ATTDEF"), scene)
if self.recenter:
self._recenter(scene, name)
elif self.pDXF is not None:
self.georeference(scene, Vector((0, 0, 0)))
if type(self.pScene) is TransverseMercator:
scene['SRID'] = "tmerc"
elif self.pScene is not None: # assume Proj
scene['SRID'] = re.findall("\+init=(.+)\s", self.pScene.srs)[0]
bpy.context.screen.scene = scene
return self.errors
# trying to import dimensions:
# self.separated_objects((block for block in self.dwg.blocks if block.name.startswith("*")))
def to_latlon(geojson, origin_espg=None):
"""
Convert a given geojson to wgs84. The original epsg must be included insde the crs
tag of geojson
"""
if isinstance(geojson, dict):
# get epsg code:
if origin_espg:
code = origin_espg
else:
code = epsg_code(geojson)
if code:
origin = Proj(init='epsg:%s' % code)
wgs84 = Proj(init='epsg:4326')
wrapped = test_wrap_coordinates(geojson['coordinates'], origin, wgs84)
new_coords = convert_coordinates(geojson['coordinates'], origin, wgs84, wrapped)
if new_coords:
geojson['coordinates'] = new_coords
try:
del geojson['crs']
except KeyError:
pass
return geojson
def __init__(self, filename):
self.filename = filename
self.file = netCDF4.Dataset(self.filename, 'r')
proj = pyproj.Proj(self.proj)
if "latitude" not in self.file.variables or "longitude" not in self.file.variables:
x, y = np.meshgrid(self.x, self.y)
self._lons, self._lats = proj(x, y, inverse=True)
else:
self._lats = self.file.variables["latitude"][:]
self._lons = self.file.variables["longitude"][:]
self._field_cache = dict()
def get_extent_projected(self, pyproj_obj):
""" Return raster extent in a projected coordinate system.
This is particularly useful for converting raster extent to the
coordinate system of a Basemap instance.
:param pyproj_obj: The projection system to convert into.
:type pyproj_obj: pyproj.Proj
:returns: extent in requested coordinate system (left, right, bottom, top)
:type: tuple
:Example:
>>> from mpl_toolkits.basemap import Basemap
>>> my_im = georaster.SingleBandRaster('myfile.tif')
>>> my_map = Basemap(...)
>>> my_map.imshow(my_im.r,extent=my_im.get_extent_basemap(my_map))
"""
if self.proj != None:
xll,xur,yll,yur = self.get_extent_latlon()
else:
xll,xur,yll,yur = self.extent
left,bottom = pyproj_obj(xll,yll)
right,top = pyproj_obj(xur,yur)
return (left,right,bottom,top)
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 georeference(self, scene, center):
if "latitude" not in scene and "longitude" not in scene:
if type(self.pScene) is TransverseMercator:
scene['latitude'] = self.pScene.lat
scene['longitude'] = self.pScene.lon
scene['altitude'] = 0
elif type(self.pScene) is not None:
wgs84 = Proj(init="EPSG:4326")
latlon = transform(self.pScene, wgs84, center[0], center[1], center[2])
scene['latitude'] = latlon[0]
scene['longitude'] = latlon[1]
scene['altitude'] = latlon[2]
def entities(self, name, scene=None):
"""
Iterates over all DXF entities according to the options set by user.
"""
if scene is None:
scene = bpy.context.scene
self.current_scene = scene
if self.recenter:
self.objects_before += scene.objects[:]
if self.combination != SEPARATED:
self.combined_objects((en for en in self.dwg.modelspace() if is_.combined_entity(en)), scene)
self.separated_entities((en for en in self.dwg.modelspace() if is_.separated_entity(en)), scene)
else:
self.separated_entities((en for en in self.dwg.modelspace() if en.dxftype != "ATTDEF"), scene)
if self.recenter:
self._recenter(scene, name)
elif self.pDXF is not None:
self.georeference(scene, Vector((0, 0, 0)))
if type(self.pScene) is TransverseMercator:
scene['SRID'] = "tmerc"
elif self.pScene is not None: # assume Proj
scene['SRID'] = re.findall("\+init=(.+)\s", self.pScene.srs)[0]
bpy.context.screen.scene = scene
return self.errors
# trying to import dimensions:
# self.separated_objects((block for block in self.dwg.blocks if block.name.startswith("*")))
def _get_extents(self, proj, res, lon_arr, lat_arr):
p = Proj(proj)
res = float(res)
first_good = self._first_good_nav(lon_arr, lat_arr)
one_x, one_y = p(lon_arr[first_good], lat_arr[first_good])
left_x = one_x - res * first_good[1]
right_x = left_x + res * lon_arr.shape[1]
top_y = one_y + res * first_good[0]
bot_y = top_y - res * lon_arr.shape[0]
half_x = res / 2.
half_y = res / 2.
return (left_x - half_x,
bot_y - half_y,
right_x + half_x,
top_y + half_y)
def test_lonlat_from_geos(self):
"""Get lonlats from geos."""
geos_area = mock.MagicMock()
lon_0 = 0
h = 35785831.00
geos_area.proj_dict = {'a': 6378169.00,
'b': 6356583.80,
'h': h,
'lon_0': lon_0}
expected = np.array((lon_0, 0))
import pyproj
proj = pyproj.Proj(proj='geos', **geos_area.proj_dict)
expected = proj(0, 0, inverse=True)
np.testing.assert_allclose(expected,
hf._lonlat_from_geos_angle(0, 0, geos_area))
expected = proj(0, 1000000, inverse=True)
np.testing.assert_allclose(expected,
hf._lonlat_from_geos_angle(0, 1000000 / h,
geos_area))
expected = proj(1000000, 0, inverse=True)
np.testing.assert_allclose(expected,
hf._lonlat_from_geos_angle(1000000 / h, 0,
geos_area))
expected = proj(2000000, -2000000, inverse=True)
np.testing.assert_allclose(expected,
hf._lonlat_from_geos_angle(2000000 / h,
-2000000 / h,
geos_area))
def _fill_sector_info(self):
for sector_info in self.scmi_sectors.values():
p = Proj(sector_info['projection'])
if 'lower_left_xy' in sector_info:
sector_info['lower_left_lonlat'] = p(*sector_info['lower_left_xy'], inverse=True)
else:
sector_info['lower_left_xy'] = p(*sector_info['lower_left_lonlat'])
if 'upper_right_xy' in sector_info:
sector_info['upper_right_lonlat'] = p(*sector_info['upper_right_xy'], inverse=True)
else:
sector_info['upper_right_xy'] = p(*sector_info['upper_right_lonlat'])
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)
