我们从Python开源项目中，提取了以下43个代码示例，用于说明如何使用pyproj.transform()。
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 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 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 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
def __getitem__(self, geometry): if isinstance(geometry, BaseGeometry) or getattr(geometry, "__geo_interface__", None) is not None: if self._tms_meta._bounds is None: return self.aoi(geojson=mapping(geometry), from_proj=self.proj) image = GeoImage.__getitem__(self, geometry) image._tms_meta = self._tms_meta return image else: result = super(TmsImage, self).__getitem__(geometry) image = super(TmsImage, 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__ image._tms_meta = self._tms_meta return image
def affine(self): """ The geo transform of the image Returns: affine (dict): The image's affine transform """ # TODO add check for Ratpoly or whatevs return self.__geo_transform__._affine
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 __contains__(self, g): geometry = ops.transform(self.__geo_transform__.rev, g) img_bounds = box(0, 0, *self.shape[2:0:-1]) return img_bounds.contains(geometry)
def get_tile_geometry(path, origin_espg, tolerance=500): """ Calculate the data and tile geometry for sentinel-2 tiles """ with rasterio.open(path) as src: # Get tile geometry b = src.bounds tile_shape = Polygon([(b[0], b[1]), (b[2], b[1]), (b[2], b[3]), (b[0], b[3]), (b[0], b[1])]) tile_geojson = mapping(tile_shape) # read first band of the image image = src.read(1) # create a mask of zero values mask = image == 0. # generate shapes of the mask novalue_shape = shapes(image, mask=mask, transform=src.affine) # generate polygons using shapely novalue_shape = [Polygon(s['coordinates'][0]) for (s, v) in novalue_shape] if novalue_shape: # Make sure polygons are united # also simplify the resulting polygon union = cascaded_union(novalue_shape) # generates a geojson data_shape = tile_shape.difference(union) # If there are multipolygons, select the largest one if data_shape.geom_type == 'MultiPolygon': areas = {p.area: i for i, p in enumerate(data_shape)} largest = max(areas.keys()) data_shape = data_shape[areas[largest]] # if the polygon has interior rings, remove them if list(data_shape.interiors): data_shape = Polygon(data_shape.exterior.coords) data_shape = data_shape.simplify(tolerance, preserve_topology=False) data_geojson = mapping(data_shape) else: data_geojson = tile_geojson # convert cooridnates to degrees return (to_latlon(tile_geojson, origin_espg), to_latlon(data_geojson, origin_espg))
def getNoSampleGrid(yespoints, xvar, yvar, dx, h1, h2): """Return the grid from which "no" pixels can successfully be sampled. Args: yespoints: Sequence of (x,y) points (meters) where landslide/liquefaction was observed. xvar: Numpy array of centers of columns of sampling grid. yvar: Numpy array of centers of rows of sampling grid. dx: Sampling resolution in x and y (meters). h1: Minimum buffer size for sampling non-hazard points. h2: Maximum buffer size for sampling non-hazard points. Returns: Grid of shape (len(yvar),len(xvar)) where 1's represent pixels from which "no" values can be sampled. """ shp = (len(xvar), len(yvar)) west = xvar.min() - dx / 2.0 # ?? north = yvar.max() + dx / 2.0 # ?? affine = affine_from_corner(west, north, dx, dx) donuts = [] holes = [] for h, k in yespoints: donut = createCirclePolygon(h, k, h2, dx) hole = createCirclePolygon(h, k, h1, dx) donuts.append(donut) holes.append(hole) donutburn = ((mapping(g), 1) for g in donuts) holeburn = ((mapping(g), 2) for g in holes) # we only want those pixels set where the polygon encloses the center point alltouched = False donutimg = rasterio.features.rasterize(donutburn, out_shape=shp, transform=affine, all_touched=alltouched) holeimg = rasterio.features.rasterize(holeburn, out_shape=shp, transform=affine, all_touched=alltouched) holeimg[holeimg == 0] = 1 holeimg[holeimg == 2] = 0 sampleimg = np.bitwise_and(donutimg, holeimg) return sampleimg
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 trilaterate(points): LatA = points[0][1] LonA = points[0][0] DistA = points[0][2]*MagicDistanceFactor LatB = points[1][1] LonB = points[1][0] DistB = points[1][2]*MagicDistanceFactor LatC = points[2][1] LonC = points[2][0] DistC = points[2][2]*MagicDistanceFactor # Transform from Latitude/Longitude to ECEF coordinates. xA,yA,zA = pyproj.transform(lla, ecef, LonA, LatA, 0, radians=False) xB,yB,zB = pyproj.transform(lla, ecef, LonB, LatB, 0, radians=False) xC,yC,zC = pyproj.transform(lla, ecef, LonC, LatC, 0, radians=False) # Convert to numpy arrays. P1 = numpy.array([xA/1000.0, yA/1000.0, zA/1000.0]) P2 = numpy.array([xB/1000.0, yB/1000.0, zB/1000.0]) P3 = numpy.array([xC/1000.0, yC/1000.0, zC/1000.0]) # Sphere intersection from Wikipedia + Stackoverflow. ex = (P2 - P1)/(numpy.linalg.norm(P2 - P1)) i = numpy.dot(ex, P3 - P1) ey = (P3 - P1 - i*ex)/(numpy.linalg.norm(P3 - P1 - i*ex)) ez = numpy.cross(ex,ey) d = numpy.linalg.norm(P2 - P1) j = numpy.dot(ey, P3 - P1) x = (pow(DistA,2) - pow(DistB,2) + pow(d,2))/(2*d) y = ((pow(DistA,2) - pow(DistC,2) + pow(i,2) + pow(j,2))/(2*j)) - ((i/j)*x) # Handle errors. if pow(DistA,2) - pow(x,2) - pow(y,2) < 0: return [] z = numpy.sqrt(pow(DistA,2) - pow(x,2) - pow(y,2)) lon,lat,altitude = pyproj.transform(ecef, lla, x*1000,y*1000,z*1000, radians=False) lon,lat,altitude = pyproj.transform(ecef, lla, x*1000,y*1000,z*1000, radians=True) triPt = P1 + x*ex + y*ey + z*ez lon,lat,altitude = pyproj.transform(ecef, lla, triPt[0]*1000,triPt[1]*1000,triPt[2]*1000, radians=False) return [lat,lon] # This was copied from trilat_optproblem.py and changed for only three point # combinations.
def reproject_dataset_example(dataset, dataset_example, method=1): # open dataset that must be transformed try: if dataset.split('.')[-1] == 'tif': g = gdal.Open(dataset) else: g = dataset except: g = dataset epsg_from = Get_epsg(g) # open dataset that is used for transforming the dataset try: if dataset_example.split('.')[-1] == 'tif': gland = gdal.Open(dataset_example) else: gland = dataset_example except: gland = dataset_example epsg_to = Get_epsg(gland) # Set the EPSG codes osng = osr.SpatialReference() osng.ImportFromEPSG(epsg_to) wgs84 = osr.SpatialReference() wgs84.ImportFromEPSG(epsg_from) # Get shape and geo transform from example geo_land = gland.GetGeoTransform() col=gland.RasterXSize rows=gland.RasterYSize # Create new raster mem_drv = gdal.GetDriverByName('MEM') dest1 = mem_drv.Create('', col, rows, 1, gdal.GDT_Float32) dest1.SetGeoTransform(geo_land) dest1.SetProjection(osng.ExportToWkt()) # Perform the projection/resampling if method is 1: gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_NearestNeighbour) if method is 2: gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_Bilinear) if method is 3: gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_Lanczos) if method is 4: gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_Average) return(dest1)
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=(',', ': '))
