Python rasterio.warp.calculate_default_transform() Examples

def ConvertRaster2LatLong(InputRasterFile,OutputRasterFile):

    """
    Convert a raster to lat long WGS1984 EPSG:4326 coordinates for global plotting

    MDH

    """

    # import modules
    import rasterio
    from rasterio.warp import reproject, calculate_default_transform as cdt, Resampling

    # read the source raster
    with rasterio.open(InputRasterFile) as src:
        #get input coordinate system
        Input_CRS = src.crs
        # define the output coordinate system
        Output_CRS = {'init': "epsg:4326"}
        # set up the transform
        Affine, Width, Height = cdt(Input_CRS,Output_CRS,src.width,src.height,*src.bounds)
        kwargs = src.meta.copy()
        kwargs.update({
            'crs': Output_CRS,
            'transform': Affine,
            'affine': Affine,
            'width': Width,
            'height': Height
        })

        with rasterio.open(OutputRasterFile, 'w', **kwargs) as dst:
            for i in range(1, src.count+1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.affine,
                    src_crs=src.crs,
                    dst_transform=Affine,
                    dst_crs=Output_CRS,
                    resampling=Resampling.bilinear) 

def get_max_zoom(self, snap=0.5, max_z=23):
        """Calculate raster max zoom level."""
        dst_affine, w, h = calculate_default_transform(
            self.crs,
            "epsg:3857",
            self.meta["width"],
            self.meta["height"],
            *self.crs_bounds
        )

        res_max = max(abs(dst_affine[0]), abs(dst_affine[4]))

        tgt_z = max_z
        mpp = 0.0

        # loop through the pyramid to file the closest z level
        for z in range(1, max_z):
            mpp = _meters_per_pixel(z, 0)

            if (mpp - ((mpp / 2) * snap)) < res_max:
                tgt_z = z
                break

        return tgt_z 

def get_max_zoom(src_dst, lat=0.0, tilesize=256):
    """
    Calculate raster max zoom level.

    Parameters
    ----------
    src: rasterio.io.DatasetReader
        Rasterio io.DatasetReader object
    lat: float, optional
        Center latitude of the dataset. This is only needed in case you want to
        apply latitude correction factor to ensure consitent maximum zoom level
        (default: 0.0).
    tilesize: int, optional
        Mercator tile size (default: 256).

    Returns
    -------
    max_zoom: int
        Max zoom level.

    """
    dst_affine, w, h = calculate_default_transform(
        src_dst.crs, "epsg:3857", src_dst.width, src_dst.height, *src_dst.bounds
    )

    native_resolution = max(abs(dst_affine[0]), abs(dst_affine[4]))

    # Correction factor for web-mercator projection latitude distortion
    latitude_correction_factor = math.cos(math.radians(lat))
    corrected_resolution = native_resolution * latitude_correction_factor

    max_zoom = zoom_for_pixelsize(corrected_resolution, tilesize=tilesize)
    return max_zoom 

def _get_vrt(src: DatasetReader, rs_method: int) -> WarpedVRT:
    from terracotta.drivers.raster_base import RasterDriver
    target_crs = RasterDriver._TARGET_CRS
    vrt_transform, vrt_width, vrt_height = calculate_default_transform(
        src.crs, target_crs, src.width, src.height, *src.bounds
    )
    vrt = WarpedVRT(
        src, crs=target_crs, resampling=rs_method, transform=vrt_transform,
        width=vrt_width, height=vrt_height
    )
    return vrt 

def get_min_zoom(self, snap=0.5, max_z=23):
        """Calculate raster min zoom level."""
        dst_affine, w, h = calculate_default_transform(
            self.crs,
            "epsg:3857",
            self.meta["width"],
            self.meta["height"],
            *self.crs_bounds
        )

        res_max = max(abs(dst_affine[0]), abs(dst_affine[4]))
        max_decim = self.overiew_levels[-1]
        resolution = max_decim * res_max

        tgt_z = 0
        mpp = 0.0

        # loop through the pyramid to file the closest z level
        for z in list(range(0, 24))[::-1]:
            mpp = _meters_per_pixel(z, 0)
            tgt_z = z

            if (mpp - ((mpp / 2) * snap)) > resolution:
                break

        return tgt_z 

def reprojectedRaster(rasterFn,ref_vectorFn):
    dst_crs=gpd.read_file(ref_vectorFn).crs
    print(dst_crs) #{'init': 'epsg:4326'}
    dst_raster_projected=os.path.join(dataFp_1,r"svf_dstRasterProjected_b.tif")
    a_T = datetime.datetime.now()
    
    # dst_crs='EPSG:4326'
    with rasterio.open(rasterFn) as src:
        transform, width, height = calculate_default_transform(src.crs, dst_crs, src.width, src.height, *src.bounds)
        kwargs = src.meta.copy()
        kwargs.update({
        'crs': dst_crs,
        'transform': transform,
        'width': width,
        'height': height,
        # 'compress': "LZW",
        'dtype':rasterio.float32,
        })
        # print(src.count)

        with rasterio.open(dst_raster_projected, 'w', **kwargs) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=transform,
                    dst_crs=dst_crs,
                    resampling=Resampling.nearest
                    )     
    
    b_T = datetime.datetime.now()
    print("reprojected time span:", b_T-a_T)    
 

#根据Polgyon统计raster栅格信息 

def reprojectedRaster(rasterFn,ref_vectorFn,dst_raster_projected):
    dst_crs=gpd.read_file(ref_vectorFn).crs
    print(dst_crs) #{'init': 'epsg:4326'}    
    a_T = datetime.datetime.now()
    
    # dst_crs='EPSG:4326'
    with rasterio.open(rasterFn) as src:
        transform, width, height = calculate_default_transform(src.crs, dst_crs, src.width, src.height, *src.bounds)
        kwargs = src.meta.copy()
        kwargs.update({
        'crs': dst_crs,
        'transform': transform,
        'width': width,
        'height': height,
        # 'compress': "LZW",
        'dtype':rasterio.uint8,  #rasterio.float32
        })
        # print(src.count)

        with rasterio.open(dst_raster_projected, 'w', **kwargs) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=transform,
                    dst_crs=dst_crs,
                    resampling=Resampling.nearest
                    )     
    
    b_T = datetime.datetime.now()
    print("reprojected time span:", b_T-a_T)    
 
#根据Polgyon统计raster栅格信息 

def get_zooms(
    src_dst: Union[DatasetReader, DatasetWriter, WarpedVRT],
    ensure_global_max_zoom: bool = False,
    tilesize: int = 256,
) -> Tuple[int, int]:
    """
    Calculate raster min/max mercator zoom level.

    Parameters
    ----------
    src_dst: rasterio.io.DatasetReader
        Rasterio io.DatasetReader object
    ensure_global_max_zoom: bool, optional
        Apply latitude correction factor to ensure max_zoom equality for global
        datasets covering different latitudes (default: False).
    tilesize: int, optional
        Mercator tile size (default: 256).

    Returns
    -------
    min_zoom, max_zoom: Tuple
        Min/Max Mercator zoom levels.

    """
    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]
    lat = center[1] if ensure_global_max_zoom else 0

    dst_affine, w, h = calculate_default_transform(
        src_dst.crs,
        constants.WEB_MERCATOR_CRS,
        src_dst.width,
        src_dst.height,
        *src_dst.bounds,
    )

    mercator_resolution = max(abs(dst_affine[0]), abs(dst_affine[4]))

    # Correction factor for web-mercator projection latitude scale change
    latitude_correction_factor = math.cos(math.radians(lat))
    adjusted_resolution = mercator_resolution * latitude_correction_factor

    max_zoom = zoom_for_pixelsize(adjusted_resolution, tilesize=tilesize)

    ovr_resolution = adjusted_resolution * max(h, w) / tilesize
    min_zoom = zoom_for_pixelsize(ovr_resolution, tilesize=tilesize)

    return (min_zoom, max_zoom) 

def get_overview_level(
    src_dst: Union[DatasetReader, DatasetWriter, WarpedVRT],
    bounds: Tuple[float, float, float, float],
    height: int,
    width: int,
    dst_crs: CRS = constants.WEB_MERCATOR_CRS,
) -> int:
    """
    Return the overview level corresponding to the tile resolution.

    Freely adapted from https://github.com/OSGeo/gdal/blob/41993f127e6e1669fbd9e944744b7c9b2bd6c400/gdal/apps/gdalwarp_lib.cpp#L2293-L2362

    Attributes
    ----------
        src_dst : rasterio.io.DatasetReader
            Rasterio io.DatasetReader object
        bounds : list
            Bounds (left, bottom, right, top) in target crs ("dst_crs").
        height : int
            Output height.
        width : int
            Output width.
        dst_crs: CRS or str, optional
            Target coordinate reference system (default "epsg:3857").

    Returns
    -------
        ovr_idx: Int or None
            Overview level

    """
    dst_transform, _, _ = calculate_default_transform(
        src_dst.crs, dst_crs, src_dst.width, src_dst.height, *src_dst.bounds
    )
    src_res = dst_transform.a

    # Compute what the "natural" output resolution
    # (in pixels) would be for this input dataset
    vrt_transform = from_bounds(*bounds, width, height)
    target_res = vrt_transform.a

    ovr_idx = -1
    if target_res > src_res:
        res = [src_res * decim for decim in src_dst.overviews(1)]

        for ovr_idx in range(ovr_idx, len(res) - 1):
            ovrRes = src_res if ovr_idx < 0 else res[ovr_idx]
            nextRes = res[ovr_idx + 1]
            if (ovrRes < target_res) and (nextRes > target_res):
                break
            if abs(ovrRes - target_res) < 1e-1:
                break
        else:
            ovr_idx = len(res) - 1

    return ovr_idx 

def get_vrt_transform(
    src_dst: Union[DatasetReader, DatasetWriter, WarpedVRT],
    bounds: Tuple[float, float, float, float],
    height: Optional[int] = None,
    width: Optional[int] = None,
    dst_crs: CRS = constants.WEB_MERCATOR_CRS,
) -> Tuple[Affine, int, int]:
    """
    Calculate VRT transform.

    Attributes
    ----------
    src_dst : rasterio.io.DatasetReader
        Rasterio io.DatasetReader object
    bounds : list
        Bounds (left, bottom, right, top) in target crs ("dst_crs").
    height : int, optional
        Desired output height of the array for the bounds.
    width : int, optional
        Desired output width of the array for the bounds.
    dst_crs: CRS or str, optional
        Target coordinate reference system (default "epsg:3857").

    Returns
    -------
    vrt_transform: Affine
        Output affine transformation matrix
    vrt_width, vrt_height: int
        Output dimensions

    """
    dst_transform, _, _ = calculate_default_transform(
        src_dst.crs, dst_crs, src_dst.width, src_dst.height, *src_dst.bounds
    )
    w, s, e, n = bounds

    if not height or not width:
        vrt_width = math.ceil((e - w) / dst_transform.a)
        vrt_height = math.ceil((s - n) / dst_transform.e)
        vrt_transform = from_bounds(w, s, e, n, vrt_width, vrt_height)
        return vrt_transform, vrt_width, vrt_height

    tile_transform = from_bounds(w, s, e, n, width, height)
    w_res = (
        tile_transform.a
        if abs(tile_transform.a) < abs(dst_transform.a)
        else dst_transform.a
    )
    h_res = (
        tile_transform.e
        if abs(tile_transform.e) < abs(dst_transform.e)
        else dst_transform.e
    )

    vrt_width = math.ceil((e - w) / w_res)
    vrt_height = math.ceil((s - n) / h_res)
    vrt_transform = from_bounds(w, s, e, n, vrt_width, vrt_height)

    return vrt_transform, vrt_width, vrt_height 

def _reproject(input_data, input_type, input_crs, target_crs, dest_path,
               resampling_method='bicubic'):

    input_crs = _check_crs(input_crs)
    target_crs = _check_crs(target_crs)
    if input_type == 'vector':
        output = input_data.to_crs(target_crs)
        if dest_path is not None:
            output.to_file(dest_path, driver='GeoJSON')

    elif input_type == 'raster':

        if isinstance(input_data, rasterio.DatasetReader):
            transform, width, height = calculate_default_transform(
                input_crs.to_wkt("WKT1_GDAL"), target_crs.to_wkt("WKT1_GDAL"),
                input_data.width, input_data.height, *input_data.bounds
            )
            kwargs = input_data.meta.copy()
            kwargs.update({'crs': target_crs.to_wkt("WKT1_GDAL"),
                           'transform': transform,
                           'width': width,
                           'height': height})

            if dest_path is not None:
                with rasterio.open(dest_path, 'w', **kwargs) as dst:
                    for band_idx in range(1, input_data.count + 1):
                        rasterio.warp.reproject(
                            source=rasterio.band(input_data, band_idx),
                            destination=rasterio.band(dst, band_idx),
                            src_transform=input_data.transform,
                            src_crs=input_data.crs,
                            dst_transform=transform,
                            dst_crs=target_crs.to_wkt("WKT1_GDAL"),
                            resampling=getattr(Resampling, resampling_method)
                        )
                output = rasterio.open(dest_path)
                input_data.close()

            else:
                output = np.zeros(shape=(height, width, input_data.count))
                for band_idx in range(1, input_data.count + 1):
                    rasterio.warp.reproject(
                        source=rasterio.band(input_data, band_idx),
                        destination=output[:, :, band_idx-1],
                        src_transform=input_data.transform,
                        src_crs=input_data.crs,
                        dst_transform=transform,
                        dst_crs=target_crs,
                        resampling=getattr(Resampling, resampling_method)
                    )

        elif isinstance(input_data, gdal.Dataset):
            if dest_path is not None:
                gdal.Warp(dest_path, input_data,
                          dstSRS='EPSG:' + str(target_crs.to_epsg()))
                output = gdal.Open(dest_path)
            else:
                raise ValueError('An output path must be provided for '
                                 'reprojecting GDAL datasets.')
    return output