Python ops.transform方法代码示例

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

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

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

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

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

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

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

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

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

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

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

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

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

# 需要导入模块: from shapely import ops [as 别名]
# 或者: from shapely.ops import transform [as 别名]
def toMeters(geometry):
    project = partial(
    pyproj.transform,
    pyproj.Proj(init='EPSG:4326'),
    pyproj.Proj(init='EPSG:32633'))
    return transform(project,geometry).length 

# 需要导入模块: from shapely import ops [as 别名]
# 或者: from shapely.ops import transform [as 别名]
def _tile_coords(self, bounds):
        """ convert mercator bbox to tile index limits """
        tfm = pyproj.Transformer.from_crs(3857, 4326, always_xy=True)
        bounds = ops.transform(tfm.transform, box(*bounds)).bounds

        # because tiles have a common corner, the tiles that cover a
        # given tile includes the adjacent neighbors.
        # https://github.com/mapbox/mercantile/issues/84#issuecomment-413113791

        west, south, east, north = bounds
        epsilon = 1.0e-10
        if east != west and north != south:
            # 2D bbox
            # shrink the bounds a small amount so that
            # shapes/tiles round trip.
            west += epsilon
            south += epsilon
            east -= epsilon
            north -= epsilon

        params = [west, south, east, north, [self.zoom_level]]
        tile_coords = [(tile.x, tile.y) for tile in mercantile.tiles(*params)]
        xtiles, ytiles = zip(*tile_coords)
        minx = min(xtiles)
        miny = min(ytiles)
        maxx = max(xtiles) 
        maxy = max(ytiles)
        return minx, miny, maxx, maxy