Python rasterio.crs.CRS.from_epsg() Examples

def feature_to_mercator(feature):
    """Normalize feature and converts coords to 3857.

    Args:
      feature: geojson feature to convert to mercator geometry.
    """
    # Ref: https://gist.github.com/dnomadb/5cbc116aacc352c7126e779c29ab7abe

    src_crs = CRS.from_epsg(4326)
    dst_crs = CRS.from_epsg(3857)

    geometry = feature["geometry"]
    if geometry["type"] == "Polygon":
        xys = (zip(*part) for part in geometry["coordinates"])
        xys = (list(zip(*transform(src_crs, dst_crs, *xy))) for xy in xys)

        yield {"coordinates": list(xys), "type": "Polygon"}

    elif geometry["type"] == "MultiPolygon":
        for component in geometry["coordinates"]:
            xys = (zip(*part) for part in component)
            xys = (list(zip(*transform(src_crs, dst_crs, *xy))) for xy in xys)

            yield {"coordinates": list(xys), "type": "Polygon"} 

def get_web_optimized_params(
    src_dst,
    tilesize=256,
    latitude_adjustment: bool = True,
    warp_resampling: str = "nearest",
    grid_crs=CRS.from_epsg(3857),
) -> Dict:
    """Return VRT parameters for a WebOptimized COG."""
    bounds = list(
        transform_bounds(
            src_dst.crs, CRS.from_epsg(4326), *src_dst.bounds, densify_pts=21
        )
    )
    center = [(bounds[0] + bounds[2]) / 2, (bounds[1] + bounds[3]) / 2]

    lat = 0 if latitude_adjustment else center[1]
    max_zoom = get_max_zoom(src_dst, lat=lat, tilesize=tilesize)

    extrema = tile_extrema(bounds, max_zoom)

    left, _, _, top = mercantile.xy_bounds(
        extrema["x"]["min"], extrema["y"]["min"], max_zoom
    )
    vrt_res = _meters_per_pixel(max_zoom, 0, tilesize=tilesize)
    vrt_transform = Affine(vrt_res, 0, left, 0, -vrt_res, top)

    vrt_width = (extrema["x"]["max"] - extrema["x"]["min"]) * tilesize
    vrt_height = (extrema["y"]["max"] - extrema["y"]["min"]) * tilesize

    return dict(
        crs=grid_crs,
        transform=vrt_transform,
        width=vrt_width,
        height=vrt_height,
        resampling=ResamplingEnums[warp_resampling],
    ) 

def test_mapchete_input(mapchete_input):
    """Mapchete process as input for other process."""
    with mapchete.open(mapchete_input.path) as mp:
        config = mp.config.params_at_zoom(5)
        input_data = config["input"]["file2"]
        assert input_data.bbox()
        assert input_data.bbox(CRS.from_epsg(3857))
        mp_input = input_data.open(next(mp.get_process_tiles(5)))
        assert not mp_input.is_empty() 

def regrid(ds, dimx, dimy, **kwargs):
    """
    Interpolate Dataset or DataArray `ds` to a new grid, using rasterio's
    reproject facility.

    See also: https://mapbox.github.io/rasterio/topics/resampling.html

    Parameters
    ----------
    ds : xr.Dataset|xr.DataArray
      N-dim data on a spatial grid
    dimx : pd.Index
      New x-coordinates in destination crs.
      dimx.name MUST refer to x-coord of ds.
    dimy : pd.Index
      New y-coordinates in destination crs.
      dimy.name MUST refer to y-coord of ds.
    **kwargs :
      Arguments passed to rio.wrap.reproject; of note:
      - resampling is one of gis.Resampling.{average,cubic,bilinear,nearest}
      - src_crs, dst_crs define the different crs (default: EPSG 4326, ie latlong)
    """
    namex = dimx.name
    namey = dimy.name

    ds = maybe_swap_spatial_dims(ds, namex, namey)

    src_transform = _as_transform(ds.indexes[namex],
                                  ds.indexes[namey])
    dst_transform = _as_transform(dimx, dimy)
    dst_shape = len(dimy), len(dimx)

    kwargs.update(dst_shape=dst_shape,
                  src_transform=src_transform,
                  dst_transform=dst_transform)
    kwargs.setdefault("src_crs", CRS.from_epsg(4326))
    kwargs.setdefault("dst_crs", CRS.from_epsg(4326))

    def _reproject(src, dst_shape, **kwargs):
        dst = np.empty(src.shape[:-2] + dst_shape, dtype=src.dtype)
        rio.warp.reproject(np.asarray(src), dst, **kwargs)
        return dst

    data_vars = ds.data_vars.values() if isinstance(ds, xr.Dataset) else (ds,)
    dtypes = {da.dtype for da in data_vars}
    assert len(dtypes) == 1, "regrid can only reproject datasets with homogeneous dtype"

    return (
        xr.apply_ufunc(_reproject, ds,
                       input_core_dims=[[namey, namex]],
                       output_core_dims=[['yout', 'xout']],
                       output_dtypes=[dtypes.pop()],
                       output_sizes={'yout': dst_shape[0], 'xout': dst_shape[1]},
                       dask='parallelized',
                       kwargs=kwargs)
        .rename({'yout': namey, 'xout': namex})
        .assign_coords(**{namey: (namey, dimy, ds.coords[namey].attrs),
                            namex: (namex, dimx, ds.coords[namex].attrs)})
        .assign_attrs(**ds.attrs)
    )