Python rasterio.warp方法代码示例

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import warp [as 别名]
def ResampleRaster(InputRasterFile,OutputRasterFile,XResolution,YResolution=None,Format="ENVI"):

    """
    Description goes here...

    MDH

    """

    # import modules
    import rasterio, affine
    from rasterio.warp import reproject, Resampling

    # read the source raster
    with rasterio.open(InputRasterFile) as src:
        Array = src.read()
        OldResolution = src.res

        #setup output resolution
        if YResolution == None:
            YResolution = XResolution
        NewResolution = (XResolution,YResolution)


        # setup the transform to change the resolution
        XResRatio = OldResolution[0]/NewResolution[0]
        YResRatio = OldResolution[1]/NewResolution[1]
        NewArray = np.empty(shape=(Array.shape[0], int(round(Array.shape[1] * XResRatio)), int(round(Array.shape[2] * YResRatio))))
        Aff = src.affine
        NewAff = affine.Affine(Aff.a/XResRatio, Aff.b, Aff.c, Aff.d, Aff.e/YResRatio, Aff.f)

        # reproject the raster
        reproject(Array, NewArray, src_transform=Aff, dst_transform=NewAff, src_crs = src.crs, dst_crs = src.crs, resample=Resampling.bilinear)

        # write results to file
        with rasterio.open(OutputRasterFile, 'w', driver=src.driver, \
                            height=NewArray.shape[1],width=NewArray.shape[2], \
                            nodata=src.nodata,dtype=str(NewArray.dtype), \
                            count=src.count,crs=src.crs,transform=NewAff) as dst:
            dst.write(NewArray) 

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

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import warp [as 别名]
def reproject_dataset(geotiff_path):
    """Project a GeoTIFF to the WGS84 coordinate reference system.
    See https://mapbox.github.io/rasterio/topics/reproject.html"""

    # We want to project the GeoTIFF coordinate reference system (crs)
    # to WGS84 (e.g. into the familiar Lat/Lon pairs). WGS84 is analogous
    # to EPSG:4326
    dst_crs = 'EPSG:4326'

    with rasterio.open(geotiff_path) as src:
        transform, width, height = rasterio.warp.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
        })

        satellite_img_name = get_file_name(geotiff_path)
        out_file_name = "{}_wgs84.tif".format(satellite_img_name)
        out_path = os.path.join(WGS84_DIR, out_file_name)
        with rasterio.open(out_path, 'w', **kwargs) as dst:
            for i in range(1, src.count + 1):
                rasterio.warp.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=rasterio.warp.Resampling.nearest)

        return rasterio.open(out_path), out_path 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import warp [as 别名]
def raster_file_xyz(raster_file):
    import rasterio
    import rasterio.warp
    import mercantile

    with rasterio.open(str(raster_file)) as src:
        raster_bounds = rasterio.warp.transform_bounds(src.crs, 'epsg:4326', *src.bounds)
    raster_center_x = (raster_bounds[0] + raster_bounds[2]) / 2
    raster_center_y = (raster_bounds[1] + raster_bounds[3]) / 2

    zoom = 14
    tile = mercantile.tile(raster_center_x, raster_center_y, zoom)
    return (tile.x, tile.y, zoom) 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import warp [as 别名]
def raster_file_xyz_lowzoom(raster_file):
    import rasterio
    import rasterio.warp
    import mercantile

    with rasterio.open(str(raster_file)) as src:
        raster_bounds = rasterio.warp.transform_bounds(src.crs, 'epsg:4326', *src.bounds)
    raster_center_x = (raster_bounds[0] + raster_bounds[2]) / 2
    raster_center_y = (raster_bounds[1] + raster_bounds[3]) / 2

    zoom = 10
    tile = mercantile.tile(raster_center_x, raster_center_y, zoom)
    return (tile.x, tile.y, zoom) 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import warp [as 别名]
def get_xyz(raster_file, zoom):
    import rasterio
    import rasterio.warp
    import mercantile

    with rasterio.open(str(raster_file)) as src:
        raster_bounds = rasterio.warp.transform_bounds(src.crs, 'epsg:4326', *src.bounds)
    raster_center_x = (raster_bounds[0] + raster_bounds[2]) / 2
    raster_center_y = (raster_bounds[1] + raster_bounds[3]) / 2

    tile = mercantile.tile(raster_center_x, raster_center_y, zoom)
    return (tile.x, tile.y, zoom) 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import warp [as 别名]
def get_bounds(geo_obj, crs=None):
    """Get the ``[left, bottom, right, top]`` bounds in any CRS.

    Arguments
    ---------
    geo_obj : a georeferenced raster or vector dataset.
    crs : int, optional
        The EPSG code (or other CRS format supported by rasterio.warp)
        for the CRS the bounds should be returned in. If not provided,
        the bounds will be returned in the same crs as `geo_obj`.

    Returns
    -------
    bounds : list
        ``[left, bottom, right, top]`` bounds in the input crs (if `crs` is
        ``None``) or in `crs` if it was provided.
    """
    input_data, input_type = _parse_geo_data(geo_obj)
    if input_type == 'vector':
        bounds = list(input_data.geometry.total_bounds)
    elif input_type == 'raster':
        if isinstance(input_data, rasterio.DatasetReader):
            bounds = list(input_data.bounds)
        elif isinstance(input_data, gdal.Dataset):
            input_gt = input_data.GetGeoTransform()
            min_x = input_gt[0]
            max_x = min_x + input_gt[1]*input_data.RasterXSize
            max_y = input_gt[3]
            min_y = max_y + input_gt[5]*input_data.RasterYSize

            bounds = [min_x, min_y, max_x, max_y]

    if crs is not None:
        crs = _check_crs(crs)
        src_crs = get_crs(input_data)
        # transform bounds to desired CRS
        bounds = transform_bounds(src_crs.to_wkt("WKT1_GDAL"),
                                  crs.to_wkt("WKT1_GDAL"), *bounds)

    return bounds 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import warp [as 别名]
def _tile_worker(tile):
    """
    For each tile, and given an open rasterio src, plus a`global_args` dictionary
    with attributes of `base_val`, `interval`, and a `writer_func`,
    warp a continous single band raster to a 512 x 512 mercator tile,
    then encode this tile into RGB.

    Parameters
    -----------
    tile: list
        [x, y, z] indices of tile

    Returns
    --------
    tile, buffer
        tuple with the input tile, and a bytearray with the data encoded into
        the format created in the `writer_func`

    """
    x, y, z = tile

    bounds = [
        c
        for i in (
            mercantile.xy(*mercantile.ul(x, y + 1, z)),
            mercantile.xy(*mercantile.ul(x + 1, y, z)),
        )
        for c in i
    ]

    toaffine = transform.from_bounds(*bounds + [512, 512])

    out = np.empty((512, 512), dtype=src.meta["dtype"])

    reproject(
        rasterio.band(src, 1),
        out,
        dst_transform=toaffine,
        dst_crs="epsg:3857",
        resampling=Resampling.bilinear,
    )

    out = data_to_rgb(out, global_args["base_val"], global_args["interval"])

    return tile, global_args["writer_func"](out, global_args["kwargs"].copy(), toaffine) 

# 需要导入模块: import rasterio [as 别名]
# 或者: from rasterio import warp [as 别名]
def test_reproject():
    from rasterio.warp import reproject
    from rasterio.enums import Resampling

    with rasterio.Env():
        # As source: a 1024 x 1024 raster centered on 0 degrees E and 0
        # degrees N, each pixel covering 15".
        rows, cols = src_shape = (1024, 1024)
        # decimal degrees per pixel
        d = 1.0 / 240

        # The following is equivalent to
        # A(d, 0, -cols*d/2, 0, -d, rows*d/2).
        src_transform = rasterio.Affine.translation(
                    -cols*d/2,
                    rows*d/2) * rasterio.Affine.scale(d, -d)
        src_crs = {'init': 'EPSG:4326'}
        source = np.ones(src_shape, np.uint8) * 255

        # Destination: a 2048 x 2048 dataset in Web Mercator (EPSG:3857)
        # with origin at 0.0, 0.0.
        dst_shape = (2048, 2048)
        dst_transform = Affine.from_gdal(
            -237481.5, 425.0, 0.0, 237536.4, 0.0, -425.0)
        dst_crs = {'init': 'EPSG:3857'}
        destination = np.zeros(dst_shape, np.uint8)

        reproject(
            source,
            destination,
            src_transform=src_transform,
            src_crs=src_crs,
            dst_transform=dst_transform,
            dst_crs=dst_crs,
            resampling=Resampling.nearest)

        # Assert that the destination is only partly filled.
        assert destination.any()
        assert not destination.all()


# Testing upsample function 

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