Python rasterio.crs() Examples

The following are 22 code examples of rasterio.crs(). You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may also want to check out all available functions/classes of the module rasterio , or try the search function .
Example #1
Source File: cogeo.py    From rio-tiler with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def spatial_info(address: str) -> Dict:
    """
    Return COGEO spatial info.

    Attributes
    ----------
        address : str or PathLike object
            A dataset path or URL. Will be opened in "r" mode.

    Returns
    -------
        out : dict.

    """
    with rasterio.open(address) as src_dst:
        minzoom, maxzoom = get_zooms(src_dst)
        bounds = transform_bounds(
            src_dst.crs, constants.WGS84_CRS, *src_dst.bounds, densify_pts=21
        )
        center = [(bounds[0] + bounds[2]) / 2, (bounds[1] + bounds[3]) / 2, minzoom]

    return dict(
        address=address, bounds=bounds, center=center, minzoom=minzoom, maxzoom=maxzoom
    ) 
Example #2
Source File: geographiclib.py    From s2p with GNU Affero General Public License v3.0 6 votes vote down vote up
def pyproj_transform(x, y, in_crs, out_crs, z=None):
    """
    Wrapper around pyproj to convert coordinates from an EPSG system to another.

    Args:
        x (scalar or array): x coordinate(s), expressed in in_crs
        y (scalar or array): y coordinate(s), expressed in in_crs
        in_crs (pyproj.crs.CRS or int): input coordinate reference system or EPSG code
        out_crs (pyproj.crs.CRS or int): output coordinate reference system or EPSG code
        z (scalar or array): z coordinate(s), expressed in in_crs

    Returns:
        scalar or array: x coordinate(s), expressed in out_crs
        scalar or array: y coordinate(s), expressed in out_crs
        scalar or array (optional if z): z coordinate(s), expressed in out_crs
    """
    transformer = pyproj.Transformer.from_crs(in_crs, out_crs, always_xy=True)
    if z is None:
        return transformer.transform(x, y)
    else:
        return transformer.transform(x, y, z) 
Example #3
Source File: geographiclib.py    From s2p with GNU Affero General Public License v3.0 6 votes vote down vote up
def pyproj_crs(projparams):
    """
    Wrapper around pyproj to return a pyproj.crs.CRS object that corresponds
    to the given parameters

    Args:
        projparams (int, str, dict): CRS parameters

    Returns:
        pyproj.crs.CRS: object that defines a CRS
    """
    if isinstance(projparams, str):
        try:
            projparams = int(projparams)
        except (ValueError, TypeError):
            pass
    return pyproj.crs.CRS(projparams) 
Example #4
Source File: geographiclib.py    From s2p with GNU Affero General Public License v3.0 6 votes vote down vote up
def rasterio_crs(projparams):
    """
    Return a rasterio.crs.CRS object that corresponds to the given parameters.
    See: https://pyproj4.github.io/pyproj/stable/crs_compatibility.html#converting-from-pyproj-crs-crs-to-rasterio-crs-crs

    Args:
        projparams (int, str, dict, pyproj.CRS): PROJ parameters

    Returns:
        rasterio.crs.CRS: object that can be used with rasterio
    """
    proj_crs = pyproj_crs(projparams)
    if LooseVersion(rasterio.__gdal_version__) < LooseVersion("3.0.0"):
        rio_crs = RioCRS.from_wkt(proj_crs.to_wkt(WktVersion.WKT1_GDAL))
    else:
        rio_crs = RioCRS.from_wkt(proj_crs.to_wkt())
    return rio_crs 
Example #5
Source File: test_utils.py    From pyresample with GNU Lesser General Public License v3.0 6 votes vote down vote up
def test_cf_load_crs_from_cf_gridmapping(self):
        from pyresample.utils.cf import _load_crs_from_cf_gridmapping

        def validate_crs_nh10km(crs):
            crs_dict = crs.to_dict()
            self.assertEqual(crs_dict['proj'], 'stere')
            self.assertEqual(crs_dict['lat_0'], 90.)

        def validate_crs_llwgs84(crs):
            crs_dict = crs.to_dict()
            self.assertEqual(crs_dict['proj'], 'longlat')
            self.assertEqual(crs_dict['ellps'], 'WGS84')

        crs = _load_crs_from_cf_gridmapping(self.nc_handles['nh10km'], 'Polar_Stereographic_Grid')
        validate_crs_nh10km(crs)
        crs = _load_crs_from_cf_gridmapping(self.nc_handles['llwgs84'], 'crs')
        validate_crs_llwgs84(crs) 
Example #6
Source File: test_utils.py    From pyresample with GNU Lesser General Public License v3.0 6 votes vote down vote up
def _prepare_cf_llwgs84():
    import xarray as xr
    nlat = 19
    nlon = 37
    ds = xr.Dataset({'temp': (('lat', 'lon'), np.ma.masked_all((nlat, nlon)), {'grid_mapping': 'crs'})},
                    coords={'lat': np.linspace(-90., +90., num=nlat),
                            'lon': np.linspace(-180., +180., num=nlon)})
    ds['lat'].attrs['units'] = 'degreesN'
    ds['lat'].attrs['standard_name'] = 'latitude'
    ds['lon'].attrs['units'] = 'degreesE'
    ds['lon'].attrs['standard_name'] = 'longitude'

    ds['crs'] = 0
    ds['crs'].attrs['grid_mapping_name'] = "latitude_longitude"
    ds['crs'].attrs['longitude_of_prime_meridian'] = 0.
    ds['crs'].attrs['semi_major_axis'] = 6378137.
    ds['crs'].attrs['inverse_flattening'] = 298.257223563

    return ds 
Example #7
Source File: test_utils.py    From pyresample with GNU Lesser General Public License v3.0 6 votes vote down vote up
def test_load_area(self):
        from pyresample import load_area
        from pyresample.utils import is_pyproj2
        ease_nh = load_area(os.path.join(os.path.dirname(__file__), 'test_files', 'areas.cfg'), 'ease_nh')
        if is_pyproj2():
            # pyproj 2.0+ adds some extra parameters
            projection = ("{'R': '6371228', 'lat_0': '90', 'lon_0': '0', "
                          "'no_defs': 'None', 'proj': 'laea', 'type': 'crs', "
                          "'units': 'm', 'x_0': '0', 'y_0': '0'}")
        else:
            projection = ("{'a': '6371228.0', 'lat_0': '90.0', "
                          "'lon_0': '0.0', 'proj': 'laea', 'units': 'm'}")
        nh_str = """Area ID: ease_nh
Description: Arctic EASE grid
Projection ID: ease_nh
Projection: {}
Number of columns: 425
Number of rows: 425
Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)""".format(projection)
        self.assertEqual(nh_str, ease_nh.__str__()) 
Example #8
Source File: cogeo.py    From rio-tiler with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def bounds(address: str) -> Dict:
    """
    Retrieve image bounds.

    Attributes
    ----------
        address : str
            file url.

    Returns
    -------
        out : dict
            dictionary with image bounds.

    """
    with rasterio.open(address) as src_dst:
        bounds = transform_bounds(
            src_dst.crs, constants.WGS84_CRS, *src_dst.bounds, densify_pts=21
        )
    return dict(address=address, bounds=bounds) 
Example #9
Source File: test_utils.py    From pyresample with GNU Lesser General Public License v3.0 5 votes vote down vote up
def test_area_parser_legacy(self):
        """Test legacy area parser."""
        from pyresample import parse_area_file
        from pyresample.utils import is_pyproj2
        ease_nh, ease_sh = parse_area_file(os.path.join(os.path.dirname(__file__), 'test_files', 'areas.cfg'),
                                           'ease_nh', 'ease_sh')

        if is_pyproj2():
            # pyproj 2.0+ adds some extra parameters
            projection = ("{'R': '6371228', 'lat_0': '90', 'lon_0': '0', "
                          "'no_defs': 'None', 'proj': 'laea', 'type': 'crs', "
                          "'units': 'm', 'x_0': '0', 'y_0': '0'}")
        else:
            projection = ("{'a': '6371228.0', 'lat_0': '90.0', "
                          "'lon_0': '0.0', 'proj': 'laea', 'units': 'm'}")
        nh_str = """Area ID: ease_nh
Description: Arctic EASE grid
Projection ID: ease_nh
Projection: {}
Number of columns: 425
Number of rows: 425
Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)""".format(projection)
        self.assertEqual(ease_nh.__str__(), nh_str)
        self.assertIsInstance(ease_nh.proj_dict['lat_0'], (int, float))

        if is_pyproj2():
            projection = ("{'R': '6371228', 'lat_0': '-90', 'lon_0': '0', "
                          "'no_defs': 'None', 'proj': 'laea', 'type': 'crs', "
                          "'units': 'm', 'x_0': '0', 'y_0': '0'}")
        else:
            projection = ("{'a': '6371228.0', 'lat_0': '-90.0', "
                          "'lon_0': '0.0', 'proj': 'laea', 'units': 'm'}")
        sh_str = """Area ID: ease_sh
Description: Antarctic EASE grid
Projection ID: ease_sh
Projection: {}
Number of columns: 425
Number of rows: 425
Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)""".format(projection)
        self.assertEqual(ease_sh.__str__(), sh_str)
        self.assertIsInstance(ease_sh.proj_dict['lat_0'], (int, float)) 
Example #10
Source File: test_utils.py    From pyresample with GNU Lesser General Public License v3.0 5 votes vote down vote up
def tmptiff(width=100, height=100, transform=None, crs=None, dtype=np.uint8):
    import rasterio
    array = np.ones((width, height)).astype(dtype)
    fname = '/vsimem/%s' % uuid.uuid4()
    with rasterio.open(fname, 'w', driver='GTiff', count=1, transform=transform,
                       width=width, height=height, crs=crs, dtype=dtype) as dst:
        dst.write(array, 1)
    return fname 
Example #11
Source File: test_utils.py    From pyresample with GNU Lesser General Public License v3.0 5 votes vote down vote up
def test_proj4_str_dict_conversion(self):
        from pyresample import utils

        proj_str = "+proj=lcc +ellps=WGS84 +lon_0=-95 +no_defs"
        proj_dict = utils.proj4.proj4_str_to_dict(proj_str)
        proj_str2 = utils.proj4.proj4_dict_to_str(proj_dict)
        proj_dict2 = utils.proj4.proj4_str_to_dict(proj_str2)
        self.assertDictEqual(proj_dict, proj_dict2)
        self.assertIsInstance(proj_dict['lon_0'], float)
        self.assertIsInstance(proj_dict2['lon_0'], float)

        # EPSG
        proj_str = '+init=EPSG:4326'
        proj_dict_exp = {'init': 'EPSG:4326'}
        proj_dict = utils.proj4.proj4_str_to_dict(proj_str)
        self.assertEqual(proj_dict, proj_dict_exp)
        self.assertEqual(utils.proj4.proj4_dict_to_str(proj_dict), proj_str)  # round-trip

        proj_str = 'EPSG:4326'
        proj_dict_exp = {'init': 'EPSG:4326'}
        proj_dict_exp2 = {'proj': 'longlat', 'datum': 'WGS84', 'no_defs': None, 'type': 'crs'}
        proj_dict = utils.proj4.proj4_str_to_dict(proj_str)
        if 'init' in proj_dict:
            # pyproj <2.0
            self.assertEqual(proj_dict, proj_dict_exp)
        else:
            # pyproj 2.0+
            self.assertEqual(proj_dict, proj_dict_exp2)
        # input != output for this style of EPSG code
        # EPSG to PROJ.4 can be lossy
        # self.assertEqual(utils._proj4.proj4_dict_to_str(proj_dict), proj_str)  # round-trip 
Example #12
Source File: test_utils.py    From pyresample with GNU Lesser General Public License v3.0 5 votes vote down vote up
def test_get_area_def_from_raster(self):
        from pyresample import utils
        from rasterio.crs import CRS
        from affine import Affine
        x_size = 791
        y_size = 718
        transform = Affine(300.0379266750948, 0.0, 101985.0,
                           0.0, -300.041782729805, 2826915.0)
        crs = CRS(init='epsg:3857')
        if utils.is_pyproj2():
            # pyproj 2.0+ expands CRS parameters
            from pyproj import CRS
            proj_dict = CRS(3857).to_dict()
        else:
            proj_dict = crs.to_dict()
        source = tmptiff(x_size, y_size, transform, crs=crs)
        area_id = 'area_id'
        proj_id = 'proj_id'
        description = 'name'
        area_def = utils.rasterio.get_area_def_from_raster(
            source, area_id=area_id, name=description, proj_id=proj_id)
        self.assertEqual(area_def.area_id, area_id)
        self.assertEqual(area_def.proj_id, proj_id)
        self.assertEqual(area_def.description, description)
        self.assertEqual(area_def.width, x_size)
        self.assertEqual(area_def.height, y_size)
        self.assertDictEqual(proj_dict, area_def.proj_dict)
        self.assertTupleEqual(area_def.area_extent, (transform.c, transform.f + transform.e * y_size,
                                                     transform.c + transform.a * x_size, transform.f)) 
Example #13
Source File: utils.py    From label-maker with MIT License 5 votes vote down vote up
def get_tile_wms(tile, imagery, folder, kwargs):
    """
    Read a WMS endpoint with query parameters corresponding to a TMS tile

    Converts the tile boundaries to the spatial/coordinate reference system
    (SRS or CRS) specified by the WMS query parameter.
    """
    # retrieve the necessary parameters from the query string
    query_dict = parse_qs(imagery.lower())
    image_format = query_dict.get('format')[0].split('/')[1]
    wms_version = query_dict.get('version')[0]
    if wms_version == '1.3.0':
        wms_srs = query_dict.get('crs')[0]
    else:
        wms_srs = query_dict.get('srs')[0]

    # find our tile bounding box
    bound = bounds(*[int(t) for t in tile.split('-')])
    xmin, ymin, xmax, ymax = transform_bounds(WGS84_CRS, CRS.from_string(wms_srs), *bound, densify_pts=21)

    # project the tile bounding box from lat/lng to WMS SRS
    bbox = (
        [ymin, xmin, ymax, xmax] if wms_version == "1.3.0" else [xmin, ymin, xmax, ymax]
    )

    # request the image with the transformed bounding box and save
    wms_url = imagery.replace('{bbox}', ','.join([str(b) for b in bbox]))
    r = requests.get(wms_url, auth=kwargs.get('http_auth'))
    tile_img = op.join(folder, '{}.{}'.format(tile, image_format))
    with open(tile_img, 'wb') as w:
        w.write(r.content)
    return tile_img 
Example #14
Source File: test_utils.py    From pyresample with GNU Lesser General Public License v3.0 5 votes vote down vote up
def test_load_cf_llwgs84(self):
        from pyresample.utils import load_cf_area

        def validate_llwgs84(adef, cfinfo, lat='lat', lon='lon'):
            self.assertEqual(adef.shape, (19, 37))
            xc = adef.projection_x_coords
            yc = adef.projection_y_coords
            self.assertEqual(xc[0], -180., msg="Wrong x axis (index 0)")
            self.assertEqual(xc[1], -180. + 10.0, msg="Wrong x axis (index 1)")
            self.assertEqual(yc[0], -90., msg="Wrong y axis (index 0)")
            self.assertEqual(yc[1], -90. + 10.0, msg="Wrong y axis (index 1)")
            self.assertEqual(cfinfo['lon'], lon)
            self.assertEqual(cf_info['lat'], lat)
            self.assertEqual(cf_info['type_of_grid_mapping'], 'latitude_longitude')
            self.assertEqual(cf_info['x']['varname'], 'lon')
            self.assertEqual(cf_info['x']['first'], -180.)
            self.assertEqual(cf_info['y']['varname'], 'lat')
            self.assertEqual(cf_info['y']['first'], -90.)

        # prepare xarray Dataset
        cf_file = _prepare_cf_llwgs84()

        # load using a variable= that is a valid grid_mapping container
        adef, cf_info = load_cf_area(cf_file, 'crs', y='lat', x='lon')
        validate_llwgs84(adef, cf_info, lat=None, lon=None)

        # load using a variable=temp
        adef, cf_info = load_cf_area(cf_file, 'temp')
        validate_llwgs84(adef, cf_info)

        # load using a variable=None
        adef, cf_info = load_cf_area(cf_file)
        validate_llwgs84(adef, cf_info) 
Example #15
Source File: __init__.py    From mapchete with MIT License 5 votes vote down vote up
def read_output_metadata(metadata_json):
    params = read_json(metadata_json)
    grid = params["pyramid"]["grid"]
    if grid["type"] == "geodetic" and grid["shape"] == [2, 1]:
        warnings.warn(
            DeprecationWarning(
                "Deprecated grid shape ordering found. "
                "Please change grid shape from [2, 1] to [1, 2] in %s."
                % metadata_json
            )
        )
        params["pyramid"]["grid"]["shape"] = [1, 2]
    if "crs" in grid and isinstance(grid["crs"], str):
        crs = CRS.from_string(grid["crs"])
        warnings.warn(
            DeprecationWarning(
                "Deprecated 'srs' found in %s: '%s'. "
                "Use WKT representation instead: %s" % (
                    metadata_json, grid["crs"], pformat(dict(wkt=crs.to_wkt()))
                )
            )
        )
        params["pyramid"]["grid"].update(srs=dict(wkt=crs.to_wkt()))
    params.update(
        pyramid=BufferedTilePyramid(
            params["pyramid"]["grid"],
            metatiling=params["pyramid"].get("metatiling", 1),
            tile_size=params["pyramid"].get("tile_size", 256),
            pixelbuffer=params["pyramid"].get("pixelbuffer", 0)
        )
    )
    return params 
Example #16
Source File: geographiclib.py    From s2p with GNU Affero General Public License v3.0 5 votes vote down vote up
def crs_bbx(ll_poly, crs=None):
    """
    Compute the UTM bounding box of a given (lon, lat) polygon.

    Args:
        ll_poly ()
        crs (pyproj.crs.CRS): pyproj CRS object. If not specified, the default CRS of the UTM
            zone for the given geography is used.

    Returns:
       4-tuple with easting min/max and northing min/max
    """
    if not crs:
        utm_zone = compute_utm_zone(*ll_poly.mean(axis=0))
        epsg = epsg_code_from_utm_zone(utm_zone)
        crs = pyproj_crs(epsg)

    # convert lon lat polygon to target CRS
    easting, northing = pyproj.transform(pyproj.Proj(init="epsg:4326"),
                                         crs, ll_poly[:, 0], ll_poly[:, 1])

    # return UTM bounding box
    east_min = min(easting)
    east_max = max(easting)
    nort_min = min(northing)
    nort_max = max(northing)
    return east_min, east_max, nort_min, nort_max 
Example #17
Source File: utils.py    From label-maker with MIT License 4 votes vote down vote up
def get_tile_tif(tile, imagery, folder, kwargs):
    """
    Read a GeoTIFF with a window corresponding to a TMS tile

    The TMS tile bounds are converted to the GeoTIFF source CRS. That bounding
    box is converted to a pixel window which is read from the GeoTIFF. For
    remote files which are internally tiled, this will take advantage of HTTP
    GET Range Requests to avoid downloading the entire file. See more info at:
    http://www.cogeo.org/in-depth.html
    """
    bound = bounds(*[int(t) for t in tile.split('-')])
    imagery_offset = kwargs.get('imagery_offset') or [0, 0]
    with rasterio.open(imagery) as src:
        x_res, y_res = src.transform[0], src.transform[4]

        # offset our imagery in the "destination pixel" space
        offset_bound = dict()
        deg_per_pix_x = (bound.west - bound.east) / 256
        deg_per_pix_y = (bound.north - bound.south) / 256

        offset_bound['west'] = bound.west + imagery_offset[0] * deg_per_pix_x
        offset_bound['east'] = bound.east + imagery_offset[0] * deg_per_pix_x
        offset_bound['north'] = bound.north + imagery_offset[1] * deg_per_pix_y
        offset_bound['south'] = bound.south + imagery_offset[1] * deg_per_pix_y

        # project tile boundaries from lat/lng to source CRS
        x, y = transform(WGS84_CRS, src.crs, [offset_bound['west']], [offset_bound['north']])
        tile_ul_proj = x[0], y[0]

        x, y = transform(WGS84_CRS, src.crs, [offset_bound['east']], [offset_bound['south']])
        tile_lr_proj = x[0], y[0]

        # get origin point from the TIF
        tif_ul_proj = (src.bounds.left, src.bounds.top)

        # use the above information to calculate the pixel indices of the window
        top = int((tile_ul_proj[1] - tif_ul_proj[1]) / y_res)
        left = int((tile_ul_proj[0] - tif_ul_proj[0]) / x_res)
        bottom = int((tile_lr_proj[1] - tif_ul_proj[1]) / y_res)
        right = int((tile_lr_proj[0] - tif_ul_proj[0]) / x_res)

        window = ((top, bottom), (left, right))

        # read the first three bands (assumed RGB) of the TIF into an array
        data = np.empty(shape=(3, 256, 256)).astype(src.profile['dtype'])
        for k in (1, 2, 3):
            src.read(k, window=window, out=data[k - 1], boundless=True)

        # save
        tile_img = op.join(folder, '{}{}'.format(tile, '.jpg'))
        img = Image.fromarray(np.moveaxis(data, 0, -1), mode='RGB')
        img.save(tile_img)

    return tile_img 
Example #18
Source File: test_io.py    From mapchete with MIT License 4 votes vote down vote up
def test_reproject_geometry(landpoly):
    """Reproject geometry."""
    with fiona.open(landpoly, "r") as src:
        for feature in src:

            # WGS84 to Spherical Mercator
            out_geom = reproject_geometry(
                shape(feature["geometry"]), CRS(src.crs),
                CRS().from_epsg(3857))
            assert out_geom.is_valid

            # WGS84 to LAEA
            out_geom = reproject_geometry(
                shape(feature["geometry"]), CRS(src.crs),
                CRS().from_epsg(3035))
            assert out_geom.is_valid

            # WGS84 to WGS84
            out_geom = reproject_geometry(
                shape(feature["geometry"]), CRS(src.crs),
                CRS().from_epsg(4326))
            assert out_geom.is_valid

    # WGS84 bounds to Spherical Mercator
    big_box = box(-180, -90, 180, 90)
    reproject_geometry(big_box, CRS().from_epsg(4326), CRS().from_epsg(3857))

    # WGS84 bounds to Spherical Mercator raising clip error
    with pytest.raises(RuntimeError):
        reproject_geometry(
            big_box, CRS().from_epsg(4326), CRS().from_epsg(3857),
            error_on_clip=True
        )
    outside_box = box(-180, 87, 180, 90)
    assert reproject_geometry(
        outside_box, CRS().from_epsg(4326), CRS().from_epsg(3857),
    ).is_valid

    # empty geometry
    assert reproject_geometry(
        Polygon(), CRS().from_epsg(4326), CRS().from_epsg(3857)).is_empty
    assert reproject_geometry(
        Polygon(), CRS().from_epsg(4326), CRS().from_epsg(4326)).is_empty

    # CRS parameter
    big_box = box(-180, -90, 180, 90)
    assert reproject_geometry(
        big_box, 4326, 3857) == reproject_geometry(
        big_box, "4326", "3857")
    with pytest.raises(TypeError):
        reproject_geometry(big_box, 1.0, 1.0) 
Example #19
Source File: cogeo.py    From rio-tiler with BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
def area(
    address: str,
    bbox: Tuple[float, float, float, float],
    dst_crs: Optional[CRS] = None,
    bounds_crs: CRS = constants.WGS84_CRS,
    max_size: int = 1024,
    **kwargs: Any,
) -> Tuple[numpy.ndarray, numpy.ndarray]:

    """
    Read value from a bbox.

    Attributes
    ----------
    address: str
        file url.
    bbox: tuple
        bounds to read (left, bottom, right, top) in "bounds_crs".
    dst_crs: CRS or str, optional
        Target coordinate reference system, default is the dataset CRS.
    bounds_crs: CRS or str, optional
        bounds coordinate reference system, default is "epsg:4326"
    max_size: int, optional
        Limit output size array, default is 1024.
    kwargs: dict, optional
        These will be passed to the 'rio_tiler.reader.part' function.

    Returns
    -------
    data : numpy ndarray
    mask: numpy array

    """
    with rasterio.open(address) as src_dst:
        if not dst_crs:
            dst_crs = src_dst.crs

        return reader.part(
            src_dst,
            bbox,
            max_size=max_size,
            bounds_crs=bounds_crs,
            dst_crs=dst_crs,
            **kwargs,
        ) 
Example #20
Source File: reader.py    From rio-tiler with BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
def tile(
    src_dst: Union[DatasetReader, DatasetWriter, WarpedVRT],
    x: int,
    y: int,
    z: int,
    tilesize: int = 256,
    **kwargs,
) -> Tuple[numpy.ndarray, numpy.ndarray]:
    """
    Read mercator tile from an image.

    Attributes
    ----------
        src_dst : rasterio.io.DatasetReader
            rasterio.io.DatasetReader object
        x : int
            Mercator tile X index.
        y : int
            Mercator tile Y index.
        z : int
            Mercator tile ZOOM level.
        tilesize : int, optional
            Output tile size. Default is 256.
        kwargs : Any, optional
            Additional options to forward to part()

    Returns
    -------
        data : numpy ndarray
        mask: numpy array

    """
    bounds = transform_bounds(
        src_dst.crs, constants.WGS84_CRS, *src_dst.bounds, densify_pts=21
    )
    if not tile_exists(bounds, z, x, y):
        raise TileOutsideBounds(f"Tile {z}/{x}/{y} is outside image bounds")

    tile_bounds = mercantile.xy_bounds(mercantile.Tile(x=x, y=y, z=z))
    return part(
        src_dst,
        tile_bounds,
        tilesize,
        tilesize,
        dst_crs=constants.WEB_MERCATOR_CRS,
        **kwargs,
    ) 
Example #21
Source File: polygon.py    From solaris with Apache License 2.0 4 votes vote down vote up
def affine_transform_gdf(gdf, affine_obj, inverse=False, geom_col="geometry",
                         precision=None):
    """Perform an affine transformation on a GeoDataFrame.

    Arguments
    ---------
    gdf : :class:`geopandas.GeoDataFrame`, :class:`pandas.DataFrame`, or `str`
        A GeoDataFrame, pandas DataFrame with a ``"geometry"`` column (or a
        different column containing geometries, identified by `geom_col` -
        note that this column will be renamed ``"geometry"`` for ease of use
        with geopandas), or the path to a saved file in .geojson or .csv
        format.
    affine_obj : list or :class:`affine.Affine`
        An affine transformation to apply to `geom` in the form of an
        ``[a, b, d, e, xoff, yoff]`` list or an :class:`affine.Affine` object.
    inverse : bool, optional
        Use this argument to perform the inverse transformation.
    geom_col : str, optional
        The column in `gdf` corresponding to the geometry. Defaults to
        ``'geometry'``.
    precision : int, optional
        Decimal precision to round the geometries to. If not provided, no
        rounding is performed.
    """
    if isinstance(gdf, str):  # assume it's a geojson
        if gdf.lower().endswith('json'):
            gdf = gpd.read_file(gdf)
        elif gdf.lower().endswith('csv'):
            gdf = pd.read_csv(gdf)
        else:
            raise ValueError(
                "The file format is incompatible with this function.")
    if 'geometry' not in gdf.columns:
        gdf = gdf.rename(columns={geom_col: 'geometry'})
    if not isinstance(gdf['geometry'][0], Polygon):
        gdf['geometry'] = gdf['geometry'].apply(shapely.wkt.loads)
    gdf["geometry"] = gdf["geometry"].apply(convert_poly_coords,
                                            affine_obj=affine_obj,
                                            inverse=inverse)
    if precision is not None:
        gdf['geometry'] = gdf['geometry'].apply(
            _reduce_geom_precision, precision=precision)

    # the CRS is no longer valid - remove it
    gdf.crs = None

    return gdf 
Example #22
Source File: polygon.py    From solaris with Apache License 2.0 4 votes vote down vote up
def georegister_px_df(df, im_path=None, affine_obj=None, crs=None,
                      geom_col='geometry', precision=None, output_path=None):
    """Convert a dataframe of geometries in pixel coordinates to a geo CRS.

    Arguments
    ---------
    df : :class:`pandas.DataFrame`
        A :class:`pandas.DataFrame` with polygons in a column named
        ``"geometry"``.
    im_path : str, optional
        A filename or :class:`rasterio.DatasetReader` object containing an
        image that has the same bounds as the pixel coordinates in `df`. If
        not provided, `affine_obj` and `crs` must both be provided.
    affine_obj : `list` or :class:`affine.Affine`, optional
        An affine transformation to apply to `geom` in the form of an
        ``[a, b, d, e, xoff, yoff]`` list or an :class:`affine.Affine` object.
        Required if not using `raster_src`.
    crs : valid CRS `str`, `int`, or :class:`rasterio.crs.CRS` instance
        The coordinate reference system for the output GeoDataFrame as an EPSG
        code integer. Required if not providing a raster image to extract the
        information from.
    geom_col : str, optional
        The column containing geometry in `df`. If not provided, defaults to
        ``"geometry"``.
    precision : int, optional
        The decimal precision for output geometries. If not provided, the
        vertex locations won't be rounded.
    output_path : str, optional
        Path to save the resulting output to. If not provided, the object
        won't be saved to disk.

    """
    if im_path is not None:
        im = _check_rasterio_im_load(im_path)
        affine_obj = im.transform
        crs = im.crs
    else:
        if not affine_obj or not crs:
            raise ValueError('If an image path is not provided, '
                             'affine_obj and crs must be.')
    crs = _check_crs(crs)
    tmp_df = affine_transform_gdf(df, affine_obj, geom_col=geom_col,
                                  precision=precision)
    result = gpd.GeoDataFrame(tmp_df, crs='epsg:' + str(crs.to_epsg()))

    if output_path is not None:
        if output_path.lower().endswith('json'):
            result.to_file(output_path, driver='GeoJSON')
        else:
            result.to_csv(output_path, index=False)

    return result