Python pyproj.transform方法代码示例

# 需要导入模块: import pyproj [as 别名]
# 或者: from pyproj import transform [as 别名]
def get_unit_multiplier_from_epsg(epsg):
    # Can't find any lightweight way to do this, but I depend heavily on pyproj right now
    # Until there's a better way, get the unit by converting (1,0) from 'unit' to 'meter'
    meter_proj = pyproj.Proj(init='epsg:'+str(epsg), preserve_units=False) # Forces meter
    unit_proj = pyproj.Proj(init='epsg:'+str(epsg), preserve_units=True) # Stays in native unit
    try:
        # Due to numerical precision and the requirements for foot vs survey foot
        # Use a larger number for the transform
        scale_value = 1.0e6
        x2, y2 = pyproj.transform(unit_proj, meter_proj, scale_value, 0.0)
        return x2/scale_value
    except:
        pass
    return 0.0

# Some of the epsgs found in lidar files don't represent projected coordinate systems
# but are instead the datum or other geographic reference systems 

# 需要导入模块: import pyproj [as 别名]
# 或者: from pyproj import transform [as 别名]
def dbcoordinates(db, pier, epsgcode='epsg:4326'):

    try:
        from pyproj import Proj, transform
    except ImportError:
        print ("dbcoordinates: You need to install pyproj to use this method")
        return (0.0 , 0.0)

    startlong = dbselect(db,'PierLong','PIERS','PierID = "A2"')
    startlat = dbselect(db,'PierLat','PIERS','PierID = "A2"')
    coordsys = dbselect(db,'PierCoordinateSystem','PIERS','PierID = "A2"')
    startlong = float(startlong[0].replace(',','.'))
    startlat = float(startlat[0].replace(',','.'))
    coordsys = coordsys[0].split(',')[1].lower().replace(' ','')

    # projection 1: GK M34
    p1 = Proj(init=coordsys)
    # projection 2: WGS 84
    p2 = Proj(init='epsg:4326')
    # transform this point to projection 2 coordinates.
    lon1, lat1 = transform(p1,p2,startlong,startlat)

    return (lon1,lat1) 

# 需要导入模块: import pyproj [as 别名]
# 或者: from pyproj import transform [as 别名]
def GetUTMEastingNorthing(EPSG_string,latitude,longitude):
    """This returns the easting and northing for a given latitude and longitide

    Args:
        ESPG_string (str): The ESPG code. 326XX is for UTM north and 327XX is for UTM south
        latitude (float): The latitude in WGS84
        longitude (float): The longitude in WGS84

    Returns:
        easting,northing The easting and northing in the UTM zone of your selection

    Author:
        Simon M Mudd
    """

    #print "Yo, getting this stuff: "+EPSG_string
    # The lat long are in epsg 4326 which is WGS84
    inProj = Proj(init='epsg:4326')
    outProj = Proj(init=EPSG_string)
    ea,no = transform(inProj,outProj,longitude,latitude)

    return ea,no 

# 需要导入模块: import pyproj [as 别名]
# 或者: from pyproj import transform [as 别名]
def _convert_grid_indices_crs(self, grid_indices, old_crs, new_crs):
        if _OLD_PYPROJ:
            x2, y2 = pyproj.transform(old_crs, new_crs, grid_indices[:,1],
                                    grid_indices[:,0])
        else:
            x2, y2 = pyproj.transform(old_crs, new_crs, grid_indices[:,1],
                                      grid_indices[:,0], errcheck=True,
                                      always_xy=True)
        yx2 = np.column_stack([y2, x2])
        return yx2

    # def _convert_outer_indices_crs(self, affine, shape, old_crs, new_crs):
    #     y1, x1 = self.grid_indices(affine=affine, shape=shape)
    #     lx, _ = pyproj.transform(old_crs, new_crs,
    #                               x1, np.repeat(y1[0], len(x1)))
    #     rx, _ = pyproj.transform(old_crs, new_crs,
    #                               x1, np.repeat(y1[-1], len(x1)))
    #     __, by = pyproj.transform(old_crs, new_crs,
    #                               np.repeat(x1[0], len(y1)), y1)
    #     __, uy = pyproj.transform(old_crs, new_crs,
    #                               np.repeat(x1[-1], len(y1)), y1)
    #     return by, uy, lx, rx 

# 需要导入模块: import pyproj [as 别名]
# 或者: from pyproj import transform [as 别名]
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) 

# 需要导入模块: import pyproj [as 别名]
# 或者: from pyproj import transform [as 别名]
def _transform_proj(x1, y1, epsg_src, epsg_dst, polar=False):
    if PYPROJ_VER == 1:
        proj_src = Proj(init='epsg:{0}'.format(epsg_src))
        _log_proj_crs(proj_src, proj_desc='src', espg=epsg_src)
        proj_dst = Proj(init='epsg:{0}'.format(epsg_dst))
        _log_proj_crs(proj_dst, proj_desc='dst', espg=epsg_dst)
        x2, y2 = transform(proj_src, proj_dst, x1, y1)
    elif PYPROJ_VER == 2:
        # With PROJ 6+ input axis ordering needs honored per projection, even
        #   though always_xy should fix it (doesn't seem to work for UPS)
        crs_src = CRS.from_epsg(epsg_src)
        _log_proj_crs(crs_src, proj_desc='src', espg=epsg_src)
        crs_dst = CRS.from_epsg(epsg_dst)
        _log_proj_crs(crs_dst, proj_desc='dst', espg=epsg_dst)
        ct = Transformer.from_crs(crs_src, crs_dst, always_xy=(not polar))
        if polar:
            y2, x2 = ct.transform(y1, x1)
        else:
            x2, y2 = ct.transform(x1, y1)
    else:
        raise MgrsException('pyproj version unsupported')

    return x2, y2 

# 需要导入模块: import pyproj [as 别名]
# 或者: from pyproj import transform [as 别名]
def get_area_acres(geometry):
    """ Calculate area in acres for a GeoJSON geometry
    :param geometry: A GeoJSON Polygon geometry
    :returns: Area in acres
    """

    shapely_geometry = shape(geometry)
    geom_aea = transform(
        partial(
            pyproj.transform,
            pyproj.Proj(init="EPSG:4326"),
            pyproj.Proj(
                proj="aea",
                lat1=shapely_geometry.bounds[1],
                lat2=shapely_geometry.bounds[3],
            ),
        ),
        shapely_geometry,
    )
    return round(geom_aea.area / 4046.8564224) 

# 需要导入模块: import pyproj [as 别名]
# 或者: from pyproj import transform [as 别名]
def __geodesic_point_buffer(self, lon, lat, km):
        """
        Based on: https://gis.stackexchange.com/questions/289044/creating-buffer-circle-x-kilometers-from-point-using-python

        Args:
            lon:
            lat:
            km:

        Returns:
            a list of coordinates with radius km and center (lat,long) in this projection
        """
        proj_wgs84 = pyproj.Proj(init='epsg:4326')
        # Azimuthal equidistant projection
        aeqd_proj = '+proj=aeqd +lat_0={lat} +lon_0={lon} +x_0=0 +y_0=0'
        project = partial(
            pyproj.transform,
            pyproj.Proj(aeqd_proj.format(lat=lat, lon=lon)),
            proj_wgs84)
        buf = Point(0, 0).buffer(km * 1000)  # distance in metres
        return transform(project, buf).exterior.coords[:] 

# 需要导入模块: import pyproj [as 别名]
# 或者: from pyproj import transform [as 别名]
def __init__(self, crs=None, epsg=None):
        if crs is not None and epsg is not None:
            raise ValueError('crs and epsg cannot both be specified')

        if epsg is not None:
            crs = _crs_from_epsg(epsg)
        if crs is not None:
            import pyproj
            crs_out = _crs_from_epsg(4326)
            proj_in = pyproj.Proj(preserve_units=True, **crs)
            proj_out = pyproj.Proj(preserve_units=True, **crs_out)
            self.transformfunc = partial(pyproj.transform, proj_in, proj_out)
        else:
            self.transformfunc = None

        self._features = [] 

# 需要导入模块: import pyproj [as 别名]
# 或者: from pyproj import transform [as 别名]
def murlo_pre_transform(self, x_list, y_list, local_grid):
        """Transforms Poggio Civitate local coordinates into the EPSG: 3003 projection."""
        coords = self.make_coordinate_list(x_list, y_list)
        epsg3003_x = []
        epsg3003_y = []
        for coord in coords:
            x = coord[0]
            y = coord[1]
            if local_grid == 'poggio-civitate':
                # transform by Taylor Oshan
                epsg3003_x.append(x * 0.999221692962 + y * 0.0447248683267 + 1695135.19719)
                epsg3003_y.append(x * -0.0439247185204 + y * 0.999281902346 + 4780651.43589)
            elif local_grid == 'vescovado-di-murlo':
                # transform by Taylor Oshan
                epsg3003_x.append(x * 0.87120992587 + y * 0.486029300286 + 1694396.08449)
                epsg3003_y.append(x * -0.487297729938 + y * 0.873675651295 + 4782618.57257)
            else:
                # what the hell?
                pass
        print('Local grid transformed: ' + local_grid)
        print('EPSG:3003 x' + str(epsg3003_x))
        print('EPSG:3003 y' + str(epsg3003_y))
        return epsg3003_x, epsg3003_y 

# 需要导入模块: import pyproj [as 别名]
# 或者: from pyproj import transform [as 别名]
def geodesic_point_buffer(lat, lon, meters, envelope=False):

    # get WGS 84 proj
    proj_wgs84 = pyproj.Proj(init='epsg:4326')

    # Azimuthal equidistant projection
    aeqd_proj = '+proj=aeqd +lat_0={lat} +lon_0={lon} +x_0=0 +y_0=0'
    project = partial(
        pyproj.transform,
        pyproj.Proj(aeqd_proj.format(lat=lat, lon=lon)),
        proj_wgs84)

    buf = Point(0, 0).buffer(meters)  # distance in metres

    if envelope is True:
        geom = Polygon(transform(project, buf).exterior.coords[:]).envelope
    else:
        geom = Polygon(transform(project, buf).exterior.coords[:])

    return geom.to_wkt() 

# 需要导入模块: import pyproj [as 别名]
# 或者: from pyproj import transform [as 别名]
def projectShapes(features, toCRS):
    import pyproj
    from functools import partial
    import fiona.crs as fcrs
    from shapely.geometry import shape, mapping
    from shapely.ops import transform as shpTrans

    project = partial(
        pyproj.transform,
        pyproj.Proj(fcrs.from_epsg(4326)),
        pyproj.Proj(toCRS))

    return list(
        {'geometry': mapping(
            shpTrans(
                project,
                shape(feat['geometry']))
        )} for feat in features
        ) 

# 需要导入模块: import pyproj [as 别名]
# 或者: from pyproj import transform [as 别名]
def pyproj_transform(x, y, in_epsg, out_epsg):
    """
    Wrapper around pyproj to convert coordinates from an EPSG system
    to another.

    Args:
        x (scalar or array): x coordinate(s) of the point(s), expressed in `in_epsg`
        y (scalar or array): y coordinate(s) of the point(s), expressed in `in_epsg`
        in_epsg (int): EPSG code of the input coordinate system
        out_epsg (int): EPSG code of the output coordinate system

    Returns:
        scalar or array: x coordinate(s) of the point(s) in the output EPSG system
        scalar or array: y coordinate(s) of the point(s) in the output EPSG system
    """
    with warnings.catch_warnings():  # noisy warnings here with pyproj>=2.4.2
        warnings.filterwarnings("ignore", category=FutureWarning)
        in_proj = pyproj.Proj(init="epsg:{}".format(in_epsg))
        out_proj = pyproj.Proj(init="epsg:{}".format(out_epsg))
    return pyproj.transform(in_proj, out_proj, x, y) 

# 需要导入模块: import pyproj [as 别名]
# 或者: from pyproj import transform [as 别名]
def isInRegion(lat, longitude, regionLat, regionLong):
    '''
    determine if the point lat, longitude is in the region bounded by a polygon with vertices regionLat, regionLong
    adapted from solution at http://gis.stackexchange.com/questions/79215/determine-if-point-is-within-an-irregular-polygon-using-python
    '''

    # WGS84 datum                                                               
    wgs84 = pyproj.Proj(init='EPSG:4326')                                       

    # Albers Equal Area Conic (aea)                                             
    nplaea = pyproj.Proj("+proj=laea +lat_0=90 +lon_0=-40 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")

    # Transform polygon and test point coordinates to northern lat AEAC projection    
    poly_x, poly_y = pyproj.transform(wgs84, nplaea, regionLong, regionLat)      
    point_x, point_y = pyproj.transform(wgs84, nplaea, [longitude], [lat])

    poly_proj = Polygon(zip(poly_x,poly_y))

    testPoint = Point(point_x[0], point_y[0])
    return testPoint.within(poly_proj) 

# 需要导入模块: import pyproj [as 别名]
# 或者: from pyproj import transform [as 别名]
def get_earth_radius(lon, lat, a, b):
    """Compute radius of the earth ellipsoid at the given longitude and latitude.

    Args:
        lon: Geodetic longitude (degrees)
        lat: Geodetic latitude (degrees)
        a: Semi-major axis of the ellipsoid (meters)
        b: Semi-minor axis of the ellipsoid (meters)

    Returns:
        Earth Radius (meters)

    """
    geocent = pyproj.Proj(proj='geocent', a=a, b=b, units='m')
    latlong = pyproj.Proj(proj='latlong', a=a, b=b, units='m')
    x, y, z = pyproj.transform(latlong, geocent, lon, lat, 0.)
    return np.sqrt(x**2 + y**2 + z**2)