How to use the rasterio.warp function in rasterio

To help you get started, we’ve selected a few rasterio examples, based on popular ways it is used in public projects.

Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.

github opendatacube / datacube-core / tests / test_warp.py View on Github external
def test_rio_resampling_conversion():
    import pytest

    R = rasterio.warp.Resampling
    assert resampling_s2rio('nearest') == R.nearest
    assert resampling_s2rio('bilinear') == R.bilinear
    assert resampling_s2rio('Bilinear') == R.bilinear
    assert resampling_s2rio('mode') == R.mode
    assert resampling_s2rio('average') == R.average

    with pytest.raises(ValueError):
        resampling_s2rio('no_such_mode')

    # check is_resampling_nn
    assert is_resampling_nn('nearest') is True
    assert is_resampling_nn('Nearest') is True
    assert is_resampling_nn('average') is False
    assert is_resampling_nn('no_such_mode') is False

    assert is_resampling_nn(R.nearest) is True
github zackw / active-geolocator / maps / baseline_from_raster.py View on Github external
self.resolution  = args.resolution
        self.fuzz        = args.fuzz
        self.fuzz_deg    = fuzz_degrees
        self.north       = north
        self.south       = south
        self.west        = west
        self.east        = east

        self.lon_spacing = lon_spacing
        self.lat_spacing = lat_spacing

        self.lon         = lon
        self.lat         = lat

        mtx = np.empty((n_lat, n_lon), dtype=np.float32)
        rasterio.warp.reproject(
            rasterio.band(args.raster, 1),
            mtx,
            src_transform = args.raster.affine,
            src_crs       = args.raster.crs,
            dst_crs       = rasterio.crs.CRS({
                'proj': 'longlat', 'ellps': 'WGS84', 'datum': 'WGS84',
                'no_defs': True}),
            dst_transform = rasterio.transform.from_bounds(
                west, south, east, north,
                n_lon, n_lat),
            dst_nodata    = 0,
            resampling    = rasterio.warp.Resampling.cubic)

        mtx[mtx < 0] = 0
        mtx  = np.log1p(mtx)
        mtx -= np.amin(mtx)
github mojodna / marblecutter / skadi.py View on Github external
def paste((data_src, src_crs, src_transform), data, bounds, resampling=Resampling.lanczos):
    """ "Reproject" src data into the correct position within a larger image"""

    ((left, right), (bottom, top)) = warp.transform(SKADI_CRS, src_crs, bounds[::2], bounds[1::2])
    bounds = (left, bottom, right, top)
    height, width = data_src.shape[1:]
    transform_dst, _, _ = warp.calculate_default_transform(
        src_crs, src_crs, width, height, *bounds)

    print('paste height: ', height)
    print('paste width: ', width)
    print('paste bounds: ', bounds)

    data_dst = np.empty(
        data.shape,
        dtype=data.dtype,
    )

    nodata = _nodata(data_dst.dtype)
github sentinel-hub / eo-learn / io / eolearn / io / geopedia.py View on Github external
def _reproject(self, eopatch, src_raster):
        """
        Reprojects the raster data from Geopedia's CRS (POP_WEB) to EOPatch's CRS.
        """
        height, width = src_raster.shape

        dst_raster = np.ones((height, width), dtype=self.raster_dtype)

        src_bbox = transform_bbox(eopatch.bbox, CRS.POP_WEB)
        src_transform = rasterio.transform.from_bounds(*src_bbox, width=width, height=height)

        dst_bbox = eopatch.bbox
        dst_transform = rasterio.transform.from_bounds(*dst_bbox, width=width, height=height)

        rasterio.warp.reproject(src_raster, dst_raster,
                                src_transform=src_transform, src_crs={'init': CRS.ogc_string(CRS.POP_WEB)},
                                src_nodata=0,
                                dst_transform=dst_transform, dst_crs={'init': CRS.ogc_string(eopatch.bbox.crs)},
                                dst_nodata=self.no_data_val)

        return dst_raster
github RemotePixel / remotepixel-api / remotepixel / remotepixel / l8_ndvi.py View on Github external
MR = float(utils.landsat_mtl_extract(meta_data, 'REFLECTANCE_MULT_BAND_4'))
    AR = float(utils.landsat_mtl_extract(meta_data, 'REFLECTANCE_ADD_BAND_4'))

    band_address = f'{landsat_address}_B4.TIF'
    with rio.open(band_address) as band:
        lon_srs, lat_srs = warp.transform('EPSG:4326', band.crs, [coord[0]], [coord[1]])
        b4 = list(band.sample([(lon_srs[0], lat_srs[0])]))[0]
        b4 = reflectance.reflectance(b4, MR, AR, E, src_nodata=0)[0]

    MR = float(utils.landsat_mtl_extract(meta_data, 'REFLECTANCE_MULT_BAND_5'))
    AR = float(utils.landsat_mtl_extract(meta_data, 'REFLECTANCE_ADD_BAND_5'))

    band_address = f'{landsat_address}_B5.TIF'
    with rio.open(band_address) as band:
        lon_srs, lat_srs = warp.transform('EPSG:4326', band.crs, [coord[0]], [coord[1]])
        b5 = list(band.sample([(lon_srs[0], lat_srs[0])]))[0]
        b5 = reflectance.reflectance(b5, MR, AR, E, src_nodata=0)[0]

    if b4 * b5 > 0:
        ratio = (b5 - b4) / (b5 + b4)
    else:
        ratio = 0

    out = {
        'ndvi': ratio,
        'date': scene_params['date'],
        'cloud': utils.landsat_mtl_extract(meta_data, 'CLOUD_COVER')
    }

    return out
github mapbox / rasterio / rasterio / rio / features.py View on Github external
ys = [round(v, precision) for v in ys]
                self._xs = xs
                self._ys = ys

                # Prepare keyword arguments for shapes().
                kwargs = {'transform': transform}
                if not with_nodata:
                    kwargs['mask'] = msk

                src_basename = os.path.basename(src.name)

                # Yield GeoJSON features.
                for i, (g, val) in enumerate(
                        rasterio.features.shapes(img, **kwargs)):
                    if projection == 'geographic':
                        g = rasterio.warp.transform_geom(
                            src.crs, 'EPSG:4326', g,
                            antimeridian_cutting=True, precision=precision)
                    xs, ys = zip(*coords(g))
                    yield {
                        'type': 'Feature',
                        'id': "{0}:{1}".format(src_basename, i),
                        'properties': {
                            'val': val, 'filename': src_basename
                        },
                        'bbox': [min(xs), min(ys), max(xs), max(ys)],
                        'geometry': g
                    }
github openearth / glofrim / glofrim / grids.py View on Github external
Returns the index values within destination grid of coordinates 'x' and 'y'.
        If a coordinate ref. system is provides, 'x' and 'y' will first be transformed to the projection of self object
        before computing index values

        Parameters
        ----------
        x : float
            x value in coordinate reference system
        y : float
            y value in coordinate reference system
        src_crs : string
            rasterio crs object or string defining the crs of the x, y coordinates
        """
        if (src_crs is not None) and (src_crs is not self.crs):
            # coordinate ref systems of src and dst are different, so a transform is needed
            x, y = rasterio.warp.transform(src_crs, self.crs, x, y)
        x, y = np.atleast_1d(x), np.atleast_1d(y)
        r, c = rasterio.transform.rowcol(self.transform, x, y, **kwargs)
        r, c = np.asarray(r).astype(int), np.asarray(c).astype(int)
        inside = self._inside(r, c)
        if flat_index:
            # calculate 1d index; NOTE invalid indices have value -1
            inds = np.ones_like(inside, dtype=int)
            inds = inds * -1 # fill with -1
            if np.any(inside):
                inds[inside] = self.ravel_multi_index(r[inside], c[inside])
        else:
            r[~inside], c[~inside] = -1, -1
            inds = (r, c) 
        return inds
github opendatacube / datacube-core / datacube / api / geo_xarray.py View on Github external
def _get_mean_longitude(dataset):
    x, y = _get_spatial_dims(dataset)
    mean_lat = float(dataset[x][0] + dataset[x][-1]) / 2.
    mean_lon = float(dataset[y][0] + dataset[y][-1]) / 2.
    bounds = {'left': mean_lon, 'right': mean_lon, 'top': mean_lat, 'bottom': mean_lat}
    left, bottom, right, top = rasterio.warp.transform_bounds(str(dataset.crs), 'EPSG:4326', **bounds)
    return left
github RemotePixel / remotepixel-api / remotepixel / remotepixel / l8_ndvi.py View on Github external
def point(scene, coord):

    scene_params = utils.landsat_parse_scene_id(scene)
    meta_data = utils.landsat_get_mtl(scene)
    landsat_address = f's3://landsat-pds/{scene_params["key"]}'

    E = float(utils.landsat_mtl_extract(meta_data, 'SUN_ELEVATION'))

    MR = float(utils.landsat_mtl_extract(meta_data, 'REFLECTANCE_MULT_BAND_4'))
    AR = float(utils.landsat_mtl_extract(meta_data, 'REFLECTANCE_ADD_BAND_4'))

    band_address = f'{landsat_address}_B4.TIF'
    with rio.open(band_address) as band:
        lon_srs, lat_srs = warp.transform('EPSG:4326', band.crs, [coord[0]], [coord[1]])
        b4 = list(band.sample([(lon_srs[0], lat_srs[0])]))[0]
        b4 = reflectance.reflectance(b4, MR, AR, E, src_nodata=0)[0]

    MR = float(utils.landsat_mtl_extract(meta_data, 'REFLECTANCE_MULT_BAND_5'))
    AR = float(utils.landsat_mtl_extract(meta_data, 'REFLECTANCE_ADD_BAND_5'))

    band_address = f'{landsat_address}_B5.TIF'
    with rio.open(band_address) as band:
        lon_srs, lat_srs = warp.transform('EPSG:4326', band.crs, [coord[0]], [coord[1]])
        b5 = list(band.sample([(lon_srs[0], lat_srs[0])]))[0]
        b5 = reflectance.reflectance(b5, MR, AR, E, src_nodata=0)[0]

    if b4 * b5 > 0:
        ratio = (b5 - b4) / (b5 + b4)
    else:
        ratio = 0
github mojodna / marblecutter / marblecutter / catalogs / postgis.py View on Github external
END geom,
              CASE WHEN ST_IsEmpty(mask)
                THEN 'null'
                ELSE coalesce(ST_AsGeoJSON(mask), 'null')
              END mask
            FROM candidate_rows
        """.format(
            table=self.table,
            geometry_column=self.geometry_column,
            include_geometries=bool(include_geometries),
        )

        if bounds.crs == WGS84_CRS:
            left, bottom, right, top = bounds.bounds
        else:
            left, bottom, right, top = warp.transform_bounds(
                bounds.crs, WGS84_CRS, *bounds.bounds
            )

        connection = self._pool.getconn()
        try:
            with connection as conn, conn.cursor() as cur:
                cur.execute(
                    query,
                    {
                        "minx": left if left != Infinity else -180,
                        "miny": bottom if bottom != Infinity else -90,
                        "maxx": right if right != Infinity else 180,
                        "maxy": top if top != Infinity else 90,
                        "zoom": zoom,
                        "resolution": min(resolution),
                    },