Python rasterio.crs方法代码示例

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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