我们从Python开源项目中，提取了以下48个代码示例，用于说明如何使用pyproj.Proj()。
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 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 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 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 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)
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 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 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=(',', ': '))
