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]]
python类Proj()的实例源码
def calc_area_of_polygons(clons, clats):
"""
Calculate the area of lat-lon polygons.
We project sphere onto a flat surface using an equal area projection
and then calculate the area of flat polygon.
This is slow we should do some caching to avoid recomputing.
"""
areas = np.zeros(clons.shape[1:])
areas[:] = np.NAN
for j in range(areas.shape[0]):
for i in range(areas.shape[1]):
lats = clats[:, j, i]
lons = clons[:, j, i]
lat_centre = lats[0] + abs(lats[2] - lats[1]) / 2
lon_centre = lons[0] + abs(lons[1] - lons[0]) / 2
pa = pyproj.Proj(proj_str.format(lat_centre, lon_centre))
x, y = pa(lons, lats)
cop = {"type": "Polygon", "coordinates": [zip(x, y)]}
areas[j, i] = shape(cop).area
assert(np.sum(areas) is not np.NAN)
assert(np.min(areas) > 0)
return areas
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 LatLonToMeters(lat, lon ):
""""Converts given lat/lon in WGS84 Datum to XY in Spherical Mercator EPSG:900913
We Use Lambert Azimuthal Equal Area projection.
"""
proj = pyproj.Proj(proj='laea')
mx,my = proj(lon,lat)
return mx, my
def utm2geod (zn, datum, X, Y, Z) :
''' Convert UTM coordinates to geodetic coordinates '''
p = Proj (proj='utm', zone=zn, ellps=datum)
lon, lan = p (X, Y, inverse=True)
return lat, lon, 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
def _transpix(self, geometry, gsd, dem, proj):
xmin, ymin, xmax, ymax = geometry.bounds
x = np.linspace(xmin, xmax, num=int((xmax-xmin)/gsd))
y = np.linspace(ymax, ymin, num=int((ymax-ymin)/gsd))
xv, yv = np.meshgrid(x, y, indexing='xy')
if self.proj is None:
from_proj = "EPSG:4326"
else:
from_proj = self.proj
itfm = partial(pyproj.transform, pyproj.Proj(init=proj), pyproj.Proj(init=from_proj))
xv, yv = itfm(xv, yv) # if that works
if isinstance(dem, GeoImage):
g = box(xv.min(), yv.min(), xv.max(), yv.max())
try:
dem = dem[g].compute(get=dask.get) # read(quiet=True)
except AssertionError:
dem = 0 # guessing this is indexing by a 0 width geometry.
if isinstance(dem, np.ndarray):
dem = tf.resize(np.squeeze(dem), xv.shape, preserve_range=True, order=1, mode="edge")
return self.__geo_transform__.rev(xv, yv, z=dem, _type=np.float32)[::-1]
def _reproject(self, geometry, from_proj=None, to_proj=None):
if from_proj is None:
from_proj = self._default_proj
if to_proj is None:
to_proj = self.proj if self.proj is not None else "EPSG:4326"
tfm = partial(pyproj.transform, pyproj.Proj(init=from_proj), pyproj.Proj(init=to_proj))
return ops.transform(tfm, geometry)
def proj(self, co):
"""
:param co: coordinate
:return: transformed coordinate if self.pScene is defined
"""
if self.pScene is not None and self.pDXF is not None:
u = self.dxf_unit_scale
if len(co) == 3:
c1, c2, c3 = co
else:
c1, c2 = co
c3 = 0
if u != 1.0:
c1 *= u
c2 *= u
c3 *= u
# add
add = Vector((0, 0, 0))
if "latitude" in self.current_scene and "longitude" in self.current_scene:
if PYPROJ and type(self.pScene) not in (TransverseMercator, Indicator):
wgs84 = Proj(init="EPSG:4326")
cscn_lat = self.current_scene.get('latitude', 0)
cscn_lon = self.current_scene.get('longitude', 0)
cscn_alt = self.current_scene.get('altitude', 0)
add = Vector(transform(wgs84, self.pScene, cscn_lon, cscn_lat, cscn_alt))
# projection
newco = Vector(transform(self.pDXF, self.pScene, c1, c2, c3))
newco = newco - add
if any((c == float("inf") or c == float("-inf") for c in newco)):
self.errors.add("Projection results in +/- infinity coordinates.")
return newco
else:
u = self.dxf_unit_scale
if u != 1:
if len(co) == 3:
c1, c2, c3 = co
else:
c1, c2 = co
c3 = 0
c1 *= u
c2 *= u
c3 *= u
return Vector((c1, c2, c3))
else:
return Vector((co[0], co[1], co[2] if len(co) == 3 else 0))
def draw_pyproj(self, box, scene):
valid_dxf_srid = True
# DXF SCALE
box.prop(self, "dxf_scale")
# EPSG DXF
box.alert = (self.proj_scene != 'NONE' and (not valid_dxf_srid or self.proj_dxf == 'NONE'))
box.prop(self, "proj_dxf")
box.alert = False
if self.proj_dxf == 'USER':
try:
Proj(init=self.epsg_dxf_user)
except:
box.alert = True
valid_dxf_srid = False
box.prop(self, "epsg_dxf_user")
box.alert = False
box.separator()
# EPSG SCENE
col = box.column()
# Only info in case of pre-defined EPSG from current scene.
if self.internal_using_scene_srid:
col.enabled = False
col.prop(self, "proj_scene")
if self.proj_scene == 'USER':
try:
Proj(init=self.epsg_scene_user)
except Exception as e:
col.alert = True
col.prop(self, "epsg_scene_user")
col.alert = False
col.label("") # Placeholder.
elif self.proj_scene == 'TMERC':
col.prop(self, "merc_scene_lat", text="Lat")
col.prop(self, "merc_scene_lon", text="Lon")
else:
col.label("") # Placeholder.
col.label("") # Placeholder.
# user info
if self.proj_scene != 'NONE':
if not valid_dxf_srid:
box.label("DXF SRID not valid", icon="ERROR")
if self.proj_dxf == 'NONE':
box.label("", icon='ERROR')
box.label("DXF SRID must be set, otherwise")
if self.proj_scene == 'USER':
code = self.epsg_scene_user
else:
code = self.proj_scene
box.label('Scene SRID %r is ignored!' % code)
def plot_predictions_for_a_city(df, name_of_predictions_col, city):
'''
INPUT: DataFrame with location predictions; name of column in DataFrame that
contains the predictions; city ('City, State') to plot predictions for
OUTPUT: Bokeh map that shows the actual location of all the users predicted to
be in the selected city
'''
df_ = df[df[name_of_predictions_col] == city]
# Initialize two lists to hold all the latitudes and longitudes
all_lats = []
all_longs = []
# Pull all latitudes in 'centroid' column append to all_lats
for i in df_['centroid']:
all_lats.append(i[0])
# Pull all longitudes in 'centroid' column append to all_longs
for i in df_['centroid']:
all_longs.append(i[1])
# Initialize two lists to hold all the latitudes and longitudes
# converted to web mercator
all_x = []
all_y = []
# Convert latittudes and longitudes to web mercator x and y format
for i in xrange(len(all_lats)):
pnt = transform(
partial(
pyproj.transform,
pyproj.Proj(init='EPSG:4326'),
pyproj.Proj(init='EPSG:3857')),
Point(all_longs[i], all_lats[i]))
all_x.append(pnt.x)
all_y.append(pnt.y)
p = base_plot()
p.add_tile(STAMEN_TERRAIN)
p.circle(x=all_x, y=all_y, line_color=None, fill_color='#380474', size=15, alpha=.5)
output_file("stamen_toner_plot.html")
show(p)
def proj(self, co):
"""
:param co: coordinate
:return: transformed coordinate if self.pScene is defined
"""
if self.pScene is not None and self.pDXF is not None:
u = self.dxf_unit_scale
if len(co) == 3:
c1, c2, c3 = co
else:
c1, c2 = co
c3 = 0
if u != 1.0:
c1 *= u
c2 *= u
c3 *= u
# add
add = Vector((0, 0, 0))
if "latitude" in self.current_scene and "longitude" in self.current_scene:
if PYPROJ and type(self.pScene) not in (TransverseMercator, Indicator):
wgs84 = Proj(init="EPSG:4326")
cscn_lat = self.current_scene.get('latitude', 0)
cscn_lon = self.current_scene.get('longitude', 0)
cscn_alt = self.current_scene.get('altitude', 0)
add = Vector(transform(wgs84, self.pScene, cscn_lat, cscn_lon, cscn_alt))
# projection
newco = Vector(transform(self.pDXF, self.pScene, c1, c2, c3))
newco = newco - add
if any((c == float("inf") or c == float("-inf") for c in newco)):
self.errors.add("Projection results in +/- infinity coordinates.")
return newco
else:
u = self.dxf_unit_scale
if u != 1:
if len(co) == 3:
c1, c2, c3 = co
else:
c1, c2 = co
c3 = 0
c1 *= u
c2 *= u
c3 *= u
return Vector((c1, c2, c3))
else:
return Vector((co[0], co[1], co[2] if len(co) == 3 else 0))
def draw_pyproj(self, box, scene):
valid_dxf_srid = True
# DXF SCALE
box.prop(self, "dxf_scale")
# EPSG DXF
box.alert = (self.proj_scene != 'NONE' and (not valid_dxf_srid or self.proj_dxf == 'NONE'))
box.prop(self, "proj_dxf")
box.alert = False
if self.proj_dxf == 'USER':
try:
Proj(init=self.epsg_dxf_user)
except:
box.alert = True
valid_dxf_srid = False
box.prop(self, "epsg_dxf_user")
box.alert = False
box.separator()
# EPSG SCENE
col = box.column()
# Only info in case of pre-defined EPSG from current scene.
if self.internal_using_scene_srid:
col.enabled = False
col.prop(self, "proj_scene")
if self.proj_scene == 'USER':
try:
Proj(init=self.epsg_scene_user)
except Exception as e:
col.alert = True
col.prop(self, "epsg_scene_user")
col.alert = False
col.label("") # Placeholder.
elif self.proj_scene == 'TMERC':
col.prop(self, "merc_scene_lat", text="Lat")
col.prop(self, "merc_scene_lon", text="Lon")
else:
col.label("") # Placeholder.
col.label("") # Placeholder.
# user info
if self.proj_scene != 'NONE':
if not valid_dxf_srid:
box.label("DXF SRID not valid", icon="ERROR")
if self.proj_dxf == 'NONE':
box.label("", icon='ERROR')
box.label("DXF SRID must be set, otherwise")
if self.proj_scene == 'USER':
code = self.epsg_scene_user
else:
code = self.proj_scene
box.label('Scene SRID %r is ignored!' % code)
def project(self, lon, lat):
"""Project lon, lat into grid coordinates."""
if self.proj_name == 'LAMBERT':
ellipsoids = {
'WGS-84': 'WGS84',
'GRS-80': 'GRS80',
'WGS-72': 'WGS72',
'Australian': 'aust_SA',
'Krasovsky': 'krass',
'International': 'new_intl',
'Hayford-1909': 'intl',
'Clarke-1880': 'clrk80',
'Clarke-1866': 'clrk66',
'Airy': 'airy',
'Bessel': 'bessel',
# 'Hayford-1830':
'Sphere': 'sphere'
}
try:
ellps = ellipsoids[self.proj_ellipsoid]
except KeyError:
raise ValueError(
'Ellipsoid not supported: {}'.format(self.proj_ellipsoid))
p = Proj(proj='lcc', lat_0=self.orig_lat, lon_0=self.orig_lon,
lat_1=self.first_std_paral, lat_2=self.second_std_paral,
ellps=ellps)
elif self.proj_name == 'SIMPLE':
p = Proj(proj='eqc', lat_0=self.orig_lat, lon_0=self.orig_lon)
else:
raise ValueError('Projection not supported: {}'.format(
self.proj_name))
x, y = p(lon, lat)
x = np.array(x)
y = np.array(y)
x[np.isnan(lon)] = np.nan
y[np.isnan(lat)] = np.nan
x /= 1000.
y /= 1000.
# inverse rotation
theta = np.radians(self.map_rot)
x1 = x*np.cos(theta) + y*np.sin(theta)
y1 = -x*np.sin(theta) + y*np.cos(theta)
x = x1
y = y1
# Try to return the same type of lon, lat
if not isinstance(lon, np.ndarray):
if isinstance(lon, Iterable):
x = type(lon)(x)
else:
x = float(x)
if not isinstance(lat, np.ndarray):
if isinstance(lat, Iterable):
y = type(lat)(y)
else:
y = float(y)
return x, y
def __iter__(self):
""" Returns generator over shapefile rows.
Note:
The first column is an id field, taken from the id value of each shape
The middle values are taken from the property_schema
The last column is a string named geometry, which has the wkt value, the type is geometry_type.
"""
# These imports are nere, not at the module level, so the geo
# support can be an extra
self.start()
vfs, shp_file, layer_index = self._open_file_params()
with fiona.open(shp_file, vfs=vfs, layer=layer_index) as source:
if source.crs.get('init') != 'epsg:4326':
# Project back to WGS84
project = partial(pyproj.transform,
pyproj.Proj(source.crs, preserve_units=True),
pyproj.Proj(from_epsg('4326'))
)
else:
project = None
yield self.headers
for i,s in enumerate(source):
row_data = s['properties']
shp = asShape(s['geometry'])
row = [int(s['id'])]
for col_name, elem in row_data.items():
row.append(elem)
if project:
row.append(transform(project, shp))
else:
row.append(shp)
yield row
self.finish()
def read_records(records_csv, road_projection, record_projection, tz, col_id,
col_x, col_y, col_occurred):
"""Reads records csv, projects points, and localizes datetimes
:param records_csv: Path to CSV containing record data
:param road_projection: Projection CRS for road data
:param record_projection: Projection CRS for record data
:param tz: Timezone id for record data
:param col_id: Record id column name
:param col_x: Record x-coordinate column name
:param col_y: Record y-coordinate column name
:param col_occurred: Record occurred datetime column name
"""
# Create a function for projecting a point
project = partial(
pyproj.transform,
pyproj.Proj(record_projection),
pyproj.Proj(road_projection)
)
records = []
min_occurred = None
max_occurred = None
with open(records_csv, 'rb') as records_file:
csv_reader = csv.DictReader(records_file)
for row in csv_reader:
# Collect min and max occurred datetimes, as they'll be used later on
try:
parsed_dt = parser.parse(row[col_occurred])
# Localize datetimes that aren't timezone-aware
occurred = parsed_dt if parsed_dt.tzinfo else tz.localize(parsed_dt)
except:
# Skip the record if it has an invalid datetime
continue
if not min_occurred or occurred < min_occurred:
min_occurred = occurred
if not max_occurred or occurred > max_occurred:
max_occurred = occurred
records.append({
'id': row[col_id],
'point': transform(project, Point(float(row[col_x]), float(row[col_y]))),
'occurred': occurred
})
return records, min_occurred, max_occurred
def getGeonym(self, req, resp, query=None):
resp.status = falcon.HTTP_200
geo = None
# projections utilisées pour transformation en WGS84/Lambert93
s_srs = Proj(init='EPSG:2154')
t_srs = Proj(init='EPSG:4326')
if 'x' in req.params and 'y' in req.params:
lon,lat = transform(s_srs,t_srs,req.params['x'],req.params['y'])
query = geonym.ll2geonym(lat, lon)
elif 'lat' in req.params and 'lon' in req.params:
query = geonym.ll2geonym(float(req.params['lat']), float(req.params['lon']))
elif 'geonym' in req.params:
query = req.params['geonym']
elif 'adresse' in req.params:
r = requests.get(URL_GEOCODER+'/search', params={"q":req.params['adresse'], "autocomplete":0, "limit":1}, timeout=1)
geo = json.loads(r.text)
geo['source']=URL_GEOCODER
query = geonym.ll2geonym(geo['features'][0]['geometry']['coordinates'][1], geo['features'][0]['geometry']['coordinates'][0])
if query is not None and geonym.checkGeonym(query):
rev = None
data = geonym.geonym2ll(query)
if 'reverse' in req.params and req.params['reverse']=='yes':
r = requests.get(URL_GEOCODER+'/reverse', params={"lat":data['lat'],"lon":data['lon'],"limit":1}, timeout=1)
if r.status_code == 200:
rev = json.loads(r.text)
rev['source']=URL_GEOCODER
x,y = transform(t_srs,s_srs,data['lon'],data['lat'])
# on ne retourne les coordonnées Lambert que si on est en zone Lambert93
if y > -357823.2365 and x > 6037008.6939 and y < 1313632.3628 and x< 7230727.3772:
data['x'] = int(x)
data['y'] = int(y)
data['checksum'] = geonym.checksum(query)
geojson = {"type":"Feature",
"properties":data,
"link": "http://www.geonym.fr/visu/?g=%s" % (geonym.cleanGeonym(query),),
"params":geonym.getParams(),
"geometry":{"type":"Polygon","coordinates":[[[data['west'],data['south']],[data['east'],data['south']],[data['east'],data['north']],[data['west'],data['north']],[data['west'],data['south']]]]}}
if rev is not None:
geojson['reverse'] = rev
if geo is not None:
geojson['geocode'] = geo
resp.body = json.dumps(geojson, sort_keys=True, indent=4, separators=(',', ': '))
resp.set_header('Content-type','application/json')
else:
geojson = {
"type": "Feature",
"link": "https://github.com/geonym/geonymapi",
"params": geonym.getParams()
}
resp.status = falcon.HTTP_400
resp.set_header('Content-type', 'application/json')
resp.body = json.dumps(
geojson, sort_keys=True, indent=4, separators=(',', ': '))
