Python pyproj.transform() Examples

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[:] 

def ResampleRaster(InputRasterFile,OutputRasterFile,XResolution,YResolution=None,Format="ENVI"):

    """
    Description goes here...

    MDH

    """

    # import modules
    import rasterio, affine
    from rasterio.warp import reproject, Resampling

    # read the source raster
    with rasterio.open(InputRasterFile) as src:
        Array = src.read()
        OldResolution = src.res

        #setup output resolution
        if YResolution == None:
            YResolution = XResolution
        NewResolution = (XResolution,YResolution)


        # setup the transform to change the resolution
        XResRatio = OldResolution[0]/NewResolution[0]
        YResRatio = OldResolution[1]/NewResolution[1]
        NewArray = np.empty(shape=(Array.shape[0], int(round(Array.shape[1] * XResRatio)), int(round(Array.shape[2] * YResRatio))))
        Aff = src.affine
        NewAff = affine.Affine(Aff.a/XResRatio, Aff.b, Aff.c, Aff.d, Aff.e/YResRatio, Aff.f)

        # reproject the raster
        reproject(Array, NewArray, src_transform=Aff, dst_transform=NewAff, src_crs = src.crs, dst_crs = src.crs, resample=Resampling.bilinear)

        # write results to file
        with rasterio.open(OutputRasterFile, 'w', driver=src.driver, \
                            height=NewArray.shape[1],width=NewArray.shape[2], \
                            nodata=src.nodata,dtype=str(NewArray.dtype), \
                            count=src.count,crs=src.crs,transform=NewAff) as dst:
            dst.write(NewArray) 

def ConvertRaster2LatLong(InputRasterFile,OutputRasterFile):

    """
    Convert a raster to lat long WGS1984 EPSG:4326 coordinates for global plotting

    MDH

    """

    # import modules
    import rasterio
    from rasterio.warp import reproject, calculate_default_transform as cdt, Resampling

    # read the source raster
    with rasterio.open(InputRasterFile) as src:
        #get input coordinate system
        Input_CRS = src.crs
        # define the output coordinate system
        Output_CRS = {'init': "epsg:4326"}
        # set up the transform
        Affine, Width, Height = cdt(Input_CRS,Output_CRS,src.width,src.height,*src.bounds)
        kwargs = src.meta.copy()
        kwargs.update({
            'crs': Output_CRS,
            'transform': Affine,
            'affine': Affine,
            'width': Width,
            'height': Height
        })

        with rasterio.open(OutputRasterFile, 'w', **kwargs) as dst:
            for i in range(1, src.count+1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.affine,
                    src_crs=src.crs,
                    dst_transform=Affine,
                    dst_crs=Output_CRS,
                    resampling=Resampling.bilinear) 

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() 

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
        ) 

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) 

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 = [] 

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 

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 

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 _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 

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) 

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) 

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) 

def get_satpos(self):
        """Get actual satellite position in geodetic coordinates (WGS-84).

        Returns: Longitude [deg east], Latitude [deg north] and Altitude [m]
        """
        if self.satpos is None:
            logger.debug("Computing actual satellite position")

            try:
                # Get satellite position in cartesian coordinates
                x, y, z = self._get_satpos_cart()

                # Transform to geodetic coordinates
                geocent = pyproj.Proj(proj='geocent')
                a, b = self.get_earth_radii()
                latlong = pyproj.Proj(proj='latlong', a=a, b=b, units='m')
                lon, lat, alt = pyproj.transform(geocent, latlong, x, y, z)
            except NoValidOrbitParams as err:
                logger.warning(err)
                lon = lat = alt = None

            # Cache results
            self.satpos = lon, lat, alt

        return self.satpos 

def get_orbital_parameters(self):
        """Collect orbital parameters for this file."""
        a = float(self.nc.attrs['earth_equatorial_radius'])
        b = float(self.nc.attrs['earth_polar_radius'])
        # nominal_satellite_height seems to be from the center of the earth
        h = float(self.nc.attrs['nominal_satellite_height']) - a
        lon_0 = self.nc.attrs['sub_longitude'] * 180 / np.pi  # it's in radians?
        sc_position = self.nc['sc_position'].attrs['sc_position_center_pixel']

        # convert ECEF coordinates to lon, lat, alt
        ecef = pyproj.Proj(proj='geocent', a=a, b=b)
        lla = pyproj.Proj(proj='latlong', a=a, b=b)
        sc_position = pyproj.transform(
            ecef, lla, sc_position[0], sc_position[1], sc_position[2])

        orbital_parameters = {
            'projection_longitude': float(lon_0),
            'projection_latitude': 0.0,
            'projection_altitude': h,
            'satellite_actual_longitude': sc_position[0],
            'satellite_actual_latitude': sc_position[1],
            'satellite_actual_altitude': sc_position[2],  # meters
        }
        return orbital_parameters 

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 

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) 

def Open_array_info(filename=''):
    """
    Opening a tiff info, for example size of array, projection and transform matrix.

    Keyword Arguments:
    filename -- 'C:/file/to/path/file.tif' or a gdal file (gdal.Open(filename))
        string that defines the input tiff file or gdal file

    """
    f = gdal.Open(r"%s" %filename)
    if f is None:
        print '%s does not exists' %filename
    else:
        geo_out = f.GetGeoTransform()
        proj = f.GetProjection()
        size_X = f.RasterXSize
        size_Y = f.RasterYSize
        f = None
    return(geo_out, proj, size_X, size_Y) 

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 

def _do_transform(src, dst, lons, lats, alt):
        """Run pyproj.transform and stack the results."""
        x, y, z = transform(src, dst, lons, lats, alt)
        return np.dstack((x, y, z)) 

def reproject_geometry(geom, inproj4, out_epsg):
    '''Reproject a wkt geometry based on EPSG code

    Args:
        geom (ogr-geom): an ogr geom objecct
        inproj4 (str): a proj4 string
        out_epsg (str): the EPSG code to which the geometry should transformed

    Returns
        geom (ogr-geometry object): the transformed geometry

    '''

    geom = ogr.CreateGeometryFromWkt(geom)
    # input SpatialReference
    spatial_ref_in = osr.SpatialReference()
    spatial_ref_in.ImportFromProj4(inproj4)

    # output SpatialReference
    spatial_ref_out = osr.SpatialReference()
    spatial_ref_out.ImportFromEPSG(int(out_epsg))

    # create the CoordinateTransformation
    coord_transform = osr.CoordinateTransformation(spatial_ref_in,
                                                   spatial_ref_out)
    try:
        geom.Transform(coord_transform)
    except:
        print(' ERROR: Not able to transform the geometry')
        sys.exit()

    return geom 

def reproject_coordinates(self, in_x_vals, in_y_vals):
        """Returns lists of reprojected x and y coordinate values. """
        out_x, out_y = pyproj.transform(self.input_crs,
                                        self.output_crs,
                                        in_x_vals,
                                        in_y_vals)
        return out_x, out_y 

def xy2lonlat(self, x, y, inverse=False):
        x = np.array(x)
        y = np.array(y)
        if inverse:
            return pyproj.transform(self.coord_proj,
                                    self.grid_proj,
                                    x, y)
        else:
            return pyproj.transform(self.grid_proj,
                                    self.coord_proj,
                                    x, y) 

def LLH_to_ECEF(self, lat, lon, alt):
        ecef, llh = pyproj.Proj(proj='geocent'), pyproj.Proj(proj='latlong')
        x, y, z = pyproj.transform(llh, ecef, lon, lat, alt, radians=False)
        return x/1000000, y/1000000, z/1000000 

def LLH_to_ECEF(self, lat, lon, alt):
        ecef, llh = pyproj.Proj(proj='geocent'), pyproj.Proj(proj='latlong')
        x, y, z = pyproj.transform(llh, ecef, lon, lat, alt, radians=False)
        return x/1000000, y/1000000, z/1000000 

def toMeters(geometry):
    project = partial(
    pyproj.transform,
    pyproj.Proj(init='EPSG:4326'),
    pyproj.Proj(init='EPSG:32633'))
    return transform(project,geometry).length 

def toMeters(geometry):
    project = partial(
    pyproj.transform,
    pyproj.Proj(init='EPSG:4326'),
    pyproj.Proj(init='EPSG:32633'))
    return transform(project,geometry).length 

def toMeters(geometry):
    project = partial(
    pyproj.transform,
    pyproj.Proj(init='EPSG:4326'),
    pyproj.Proj(init='EPSG:32633'))
    return transform(project,geometry).length 

def _update_grid_centre(self):
        x, y = pyproj.transform(self.coord_proj, self.grid_proj,
                                self.longitude, self.latitude)

        self.grid_centre = [x, y, self.elevation - ((self.cell_count[2] - 1)
                                                    * self.cell_size[2]) / 2]