Python shapely.ops.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 _get_centerline_lonlat(gdir):
    """Quick n dirty solution to write the centerlines as a shapefile"""

    cls = gdir.read_pickle('centerlines')
    olist = []
    for j, cl in enumerate(cls[::-1]):
        mm = 1 if j == 0 else 0
        gs = dict()
        gs['RGIID'] = gdir.rgi_id
        gs['LE_SEGMENT'] = np.rint(np.max(cl.dis_on_line) * gdir.grid.dx)
        gs['MAIN'] = mm
        tra_func = partial(gdir.grid.ij_to_crs, crs=wgs84)
        gs['geometry'] = shp_trafo(tra_func, cl.line)
        olist.append(gs)

    return olist 

def _poly_area_approximation(way, nodes):
    """ Compute the approximated area of an irregular OSM polygon.
        see: https://arachnoid.com/area_irregular_polygon/
             https://gist.github.com/robinkraft/c6de2f988c9d3f01af3c
    """
    points = []
    for node in way['nd']:
        points.append([float(nodes[node['ref']]['lat']), float(nodes[node['ref']]['lon'])])

    approx = MultiPoint(points).convex_hull
    # http://openstreetmapdata.com/info/projections
    proj = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'),
                   pyproj.Proj(init='epsg:3857'))

    converted_approximation = transform(proj, approx)

    return converted_approximation.area 

def _poly_area_shapely(way, nodes):
    """ Compute the area of an irregular OSM polygon.
        see: https://arachnoid.com/area_irregular_polygon/
             https://gist.github.com/robinkraft/c6de2f988c9d3f01af3c
    """
    points = []
    for node in way['nd']:
        points.append([float(nodes[node['ref']]['lat']), float(nodes[node['ref']]['lon'])])

    geom = {'type': 'Polygon',
            'coordinates': [points]}

    s = shape(geom)
    # http://openstreetmapdata.com/info/projections
    proj = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'),
                   pyproj.Proj(init='epsg:3857'))

    newshape = transform(proj, s)

    return newshape.area 

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 get_grouped(self, lat_attr, lon_attr, admin, attr, agg_func):
        """
        Get aggregation value for points grouped by regions.
        Returns:
            Series of aggregated values
        """
        if attr is not None:
            data = self.data.get_column_view(attr)[0]
        else:
            data = np.ones(len(self.data))

        ids, _, _ = self.get_regions(lat_attr, lon_attr, admin)
        result = pd.Series(data, dtype=float)\
            .groupby(ids)\
            .agg(AGG_FUNCS[agg_func].transform)

        return result 

def _repr_html_(self):
        title = f"<b>{self.name} [{self.designator}] ({self.type})</b>"
        shapes = ""
        title += "<ul>"

        bounds = self.bounds
        projection = pyproj.Proj(
            proj="aea",  # equivalent projection
            lat_1=bounds[1],
            lat_2=bounds[3],
            lat_0=(bounds[1] + bounds[3]) / 2,
            lon_0=(bounds[0] + bounds[2]) / 2,
        )

        for polygon in self:
            transformer = pyproj.Transformer.from_proj(
                pyproj.Proj("epsg:4326"), projection, always_xy=True
            )
            projected_shape = transform(transformer.transform, polygon.polygon)
            title += f"<li>{polygon.lower}, {polygon.upper}</li>"
            shapes += projected_shape.simplify(1e3)._repr_svg_()
        title += "</ul>"
        no_wrap_div = '<div style="white-space: nowrap; width: 12%">{}</div>'
        return title + no_wrap_div.format(shapes) 

def window_at(self, geom, window_shape):
        """Return a subsetted window of a given size, centered on a geometry object

        Useful for generating training sets from vector training data
        Will throw a ValueError if the window is not within the image bounds

        Args:
            geom (shapely,geometry): Geometry to center the image on
            window_shape (tuple): The desired shape of the image as (height, width) in pixels.

        Returns:
            image: image object of same type
        """
        # Centroids of the input geometry may not be centered on the object.
        # For a covering image we use the bounds instead.
        # This is also a workaround for issue 387.
        y_size, x_size = window_shape[0], window_shape[1]
        bounds = box(*geom.bounds)
        px = ops.transform(self.__geo_transform__.rev, bounds).centroid
        miny, maxy = int(px.y - y_size/2), int(px.y + y_size/2)
        minx, maxx = int(px.x - x_size/2), int(px.x + x_size/2)
        _, y_max, x_max = self.shape
        if minx < 0 or miny < 0 or maxx > x_max or maxy > y_max:
            raise ValueError("Input geometry resulted in a window outside of the image")
        return self[:, miny:maxy, minx:maxx] 

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 = GeoDaskImage.__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)
            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 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 compute_xy(
        self: T, projection: Union[pyproj.Proj, crs.Projection, None] = None
    ) -> T:
        """Enrich the structure with new x and y columns computed through a
        projection of the latitude and longitude columns.

        The source projection is WGS84 (EPSG 4326).
        The default destination projection is a Lambert Conformal Conical
        projection centered on the data inside the dataframe.

        Other valid projections are available:

        - as ``pyproj.Proj`` objects;
        - as ``cartopy.crs.Projection`` objects.
        """

        if isinstance(projection, crs.Projection):
            projection = pyproj.Proj(projection.proj4_init)

        if projection is None:
            projection = pyproj.Proj(
                proj="lcc",
                ellps="WGS84",
                lat_1=self.data.latitude.min(),
                lat_2=self.data.latitude.max(),
                lat_0=self.data.latitude.mean(),
                lon_0=self.data.longitude.mean(),
            )

        transformer = pyproj.Transformer.from_proj(
            pyproj.Proj("epsg:4326"), projection, always_xy=True
        )
        x, y = transformer.transform(
            self.data.longitude.values, self.data.latitude.values,
        )

        return self.__class__(self.data.assign(x=x, y=y)) 

def calculate_padded_bounds(factor, bounds):
    min_x, min_y, max_x, max_y = bounds
    dx = 0.5 * (max_x - min_x) * (factor - 1.0)
    dy = 0.5 * (max_y - min_y) * (factor - 1.0)
    return geometry.box(min_x - dx, min_y - dy, max_x + dx, max_y + dy)


# function which returns its argument, used to assign to a function variable
# to use as a null transform. flake8 insists that it be named and not a
# lambda. 

def apply_to_all_coords(fn):
    return lambda shape: transform(fn, shape)


# returns a geometry which is the given bounds expanded by `factor`. that is,
# if the original shape was a 1x1 box, the new one will be `factor`x`factor`
# box, with the same centroid as the original box. 

def transform(self, g):
        result = ops.transform(self.project, g)
        result.tags = g.tags
        return result 

def _get_centerline_lonlat(gdir):
    """Quick n dirty solution to write the centerlines as a shapefile"""

    cls = gdir.read_pickle('centerlines')
    olist = []
    for j, cl in enumerate(cls[::-1]):
        mm = 1 if j == 0 else 0
        gs = gpd.GeoSeries()
        gs['RGIID'] = gdir.rgi_id
        gs['LE_SEGMENT'] = np.rint(np.max(cl.dis_on_line) * gdir.grid.dx)
        gs['MAIN'] = mm
        tra_func = partial(gdir.grid.ij_to_crs, crs=wgs84)
        gs['geometry'] = shp_trafo(tra_func, cl.line)
        olist.append(gs)

    return olist 

def generate_stop_points(chunks: List[LineString]) -> List[Point]:
    all_points = []

    # First point is the first node on the first line
    first_chunk = chunks[0].coords
    first_pt = [Point(p) for p in first_chunk][0]
    all_points.append(first_pt)

    # Then for all other points, we get from the end of each line
    for chunk in chunks:
        last_pt = [Point(p) for p in chunk.coords][-1]
        all_points.append(last_pt)

    # Now we need to convert to a Shapely object
    ap_ma = MultiPoint(all_points)

    # So we can reproject back out of equal area to 4326
    project = partial(
        pyproj.transform,
        pyproj.Proj(init='epsg:2163'),  # source coordinate system
        pyproj.Proj(init='epsg:4326'))  # destination coordinate system

    ap_ma_reproj = transform(project, ap_ma)  # apply projection

    # Final step will be to pull out all points back into a list
    return [p for p in ap_ma_reproj] 

def transform_geom(g1, src_proj, dest_proj):
    project = partial(
        pyproj.transform,
        pyproj.Proj(init=src_proj),
        pyproj.Proj(init=dest_proj))

    g2 = transform(project, g1)

    return g2 

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 _as_transform(x, y):
    lx, rx = x[[0, -1]]
    uy, ly = y[[0, -1]]

    dx = float(rx - lx)/float(len(x)-1)
    dy = float(uy - ly)/float(len(y)-1)

    return rio.transform.from_origin(lx, uy, dx, dy) 

def reproject_shapes(shapes, p1, p2):
    """
    Project a collection of `shapes` from one projection `p1` to
    another projection `p2`

    Projections can be given as strings or instances of pyproj.Proj.
    Special care is taken for the case where the final projection is
    of type rotated pole as handled by RotProj.
    """

    if p1 == p2:
        return shapes

    if isinstance(p1, RotProj):
        if p2 == 'latlong':
            reproject_points = lambda x,y: p1(x, y, inverse=True)
        else:
            raise NotImplementedError("`p1` can only be a RotProj if `p2` is latlong!")

    if isinstance(p2, RotProj):
        shapes = reproject(shapes, p1, 'latlong')
        reproject_points = p2
    else:
        reproject_points = partial(pyproj.transform, as_projection(p1), as_projection(p2))

    def _reproject_shape(shape):
        return transform(reproject_points, shape)

    if isinstance(shapes, pd.Series):
        return shapes.map(_reproject_shape)
    elif isinstance(shapes, dict):
        return OrderedDict((k, _reproject_shape(v)) for k, v in iteritems(shapes))
    else:
        return list(map(_reproject_shape, shapes)) 

def test_rasterio_glacier_masks(self):

        # The GIS was double checked externally with IDL.
        hef_file = get_demo_file('Hintereisferner_RGI5.shp')
        entity = gpd.read_file(hef_file).iloc[0]

        gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir)
        gis.define_glacier_region(gdir)

        # specifying a source will look for a DEN in a respective folder
        self.assertRaises(ValueError, gis.rasterio_glacier_mask,
                          gdir, source='SRTM')

        # this should work
        gis.rasterio_glacier_mask(gdir, source=None)

        # read dem mask
        with rasterio.open(gdir.get_filepath('glacier_mask'),
                           'r', driver='GTiff') as ds:
            profile = ds.profile
            data = ds.read(1).astype(profile['dtype'])

        # compare projections
        self.assertEqual(ds.width, gdir.grid.nx)
        self.assertEqual(ds.height, gdir.grid.ny)
        self.assertEqual(ds.transform[0], gdir.grid.dx)
        self.assertEqual(ds.transform[4], gdir.grid.dy)
        # orgin is center for gdir grid but corner for dem_mask, so shift
        self.assertAlmostEqual(ds.transform[2], gdir.grid.x0 - gdir.grid.dx/2)
        self.assertAlmostEqual(ds.transform[5], gdir.grid.y0 - gdir.grid.dy/2)

        # compare dem_mask size with RGI area
        mask_area_km2 = data.sum() * gdir.grid.dx**2 * 1e-6
        self.assertAlmostEqual(mask_area_km2, gdir.rgi_area_km2, 1)

        # how the mask is derived from the outlines it should always be larger
        self.assertTrue(mask_area_km2 > gdir.rgi_area_km2)

        # not sure if we want such a hard coded test, but this will fail if the
        # sample data changes but could also indicate changes in rasterio
        self.assertTrue(data.sum() == 3218) 

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 get_choropleth_regions(self) -> List[_ChoroplethRegion]:
        """Recalculate regions"""
        if self.attr_lat is None:
            # if we don't have locations we can't compute regions
            return []

        _, region_info, polygons = self.get_regions(self.attr_lat,
                                                    self.attr_lon,
                                                    self.admin_level)

        regions = []
        for _id in polygons:
            if isinstance(polygons[_id], MultiPolygon):
                # some regions consist of multiple polygons
                polys = list(polygons[_id].geoms)
            else:
                polys = [polygons[_id]]

            qpolys = [self.poly2qpoly(transform(self.deg2canvas, poly))
                      for poly in polys]
            regions.append(_ChoroplethRegion(id=_id, info=region_info[_id],
                                             qpolys=qpolys))

        self.choropleth_regions = sorted(regions, key=lambda cr: cr.id)
        self.get_agg_data()
        return self.choropleth_regions 

def quantize(self, shape, bounds):
        minx, miny, maxx, maxy = bounds

        def fn(x, y, z=None):
            xfac = self.extents / (maxx - minx)
            yfac = self.extents / (maxy - miny)
            x = xfac * (x - minx)
            y = yfac * (y - miny)
            return self._round(x), self._round(y)

        return transform(fn, shape) 

def effective_data(self):
        return self.data.transform(Domain(self.effective_variables))

    # Input 

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

def project_shape(
        self, projection: Union[pyproj.Proj, crs.Projection, None] = None
    ) -> base.BaseGeometry:
        """Returns a projected representation of the shape.

        By default, an equivalent projection is applied. Equivalent projections
        locally respect areas, which is convenient for the area attribute.

        Other valid projections are available:

        - as ``pyproj.Proj`` objects;
        - as ``cartopy.crs.Projection`` objects.

        """

        if self.shape is None:
            return None

        if isinstance(projection, crs.Projection):
            projection = pyproj.Proj(projection.proj4_init)

        if projection is None:
            bounds = self.bounds
            projection = pyproj.Proj(
                proj="aea",  # equivalent projection
                lat_1=bounds[1],
                lat_2=bounds[3],
                lat_0=(bounds[1] + bounds[3]) / 2,
                lon_0=(bounds[0] + bounds[2]) / 2,
            )

        transformer = pyproj.Transformer.from_proj(
            pyproj.Proj("epsg:4326"), projection, always_xy=True
        )
        projected_shape = transform(transformer.transform, self.shape,)

        if not projected_shape.is_valid:
            warnings.warn("The chosen projection is invalid for current shape")
        return projected_shape 

def plot(
        self, ax: GeoAxesSubplot, **kwargs
    ) -> List[Artist]:  # coverage: ignore
        """Plots the trajectory on a Matplotlib axis.

        The Flight supports Cartopy axis as well with automatic projection. If
        no projection is provided, a default `PlateCarree
        <https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html#platecarree>`_
        is applied.

        Example usage:

        .. code:: python

            from traffic.drawing import Mercator
            fig, ax = plt.subplots(1, subplot_kw=dict(projection=Mercator())
            flight.plot(ax, alpha=.5)

        .. note::
            See also `geoencode() <#traffic.core.Flight.geoencode>`_ for the
            altair equivalent.

        """

        if "projection" in ax.__dict__ and "transform" not in kwargs:
            kwargs["transform"] = PlateCarree()
        if self.shape is not None:
            return ax.plot(*self.shape.xy, **kwargs)
        return [] 

def annotate(
        self, ax: GeoAxesSubplot, **kwargs
    ) -> None:  # coverage: ignore
        if "projection" in ax.__dict__:
            kwargs["transform"] = PlateCarree()
        if "s" not in kwargs:
            kwargs["s"] = self.name
        ax.text(*np.array(self.centroid), **kwargs) 

def quantize(self, shape, bounds):
        minx, miny, maxx, maxy = bounds

        def fn(x, y, z=None):
            xfac = self.extents / (maxx - minx)
            yfac = self.extents / (maxy - miny)
            x = xfac * (x - minx)
            y = yfac * (y - miny)
            return self._round(x), self._round(y)

        return transform(fn, shape)