Python gdal.Warp方法代码示例

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import Warp [as 别名]
def scale_with_gdalwarp(array, prj_in, prj_out, dims_out, gt_in, gt_out, resample_alg):

    in_src = save_img(array, gt_in, prj_in, 'MEM', noDataValue=np.nan, fieldNames=[])
    # Get template projection, extent and resolution
    extent = [gt_out[0], gt_out[3]+gt_out[5]*dims_out[0],
              gt_out[0]+gt_out[1]*dims_out[1], gt_out[3]]

    # Resample with GDAL warp
    outDs = gdal.Warp('',
                      in_src,
                      dstSRS=prj_out,
                      xRes=gt_out[1],
                      yRes=gt_out[5],
                      outputBounds=extent,
                      resampleAlg=resample_alg,
                      format="MEM")

    array_out = outDs.GetRasterBand(1).ReadAsArray()

    return array_out 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import Warp [as 别名]
def crop2grid(self, delete_temp=True):
        """
        Crops raster to extent of the geomodel grid.
        """
        cornerpoints_geo = self._get_cornerpoints(self.regular_grid_extent)
        cornerpoints_dtm = self._get_cornerpoints(self.extent)

        # self.check()

        if np.any(cornerpoints_geo[:2] - cornerpoints_dtm[:2]) != 0:
            path_dest = '_cropped_DEM.tif'
            new_bounds = (self.regular_grid_extent[[0, 2, 1, 3]])
            gdal.Warp(path_dest, self.dem, options=gdal.WarpOptions(
                options=['outputBounds'], outputBounds=new_bounds))

            self.dem = gdal.Open(path_dest)
            self.dem_zval = self.dem.ReadAsArray()
            self._get_raster_dimensions()

        print('Cropped raster to geo_model.grid.extent.') 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import Warp [as 别名]
def resample(self, new_xres, new_yres, save_path):
        """
        Decrease the pixel size of the raster.

        Args:
            new_xres (int): desired resolution in x-direction
            new_yres (int): desired resolution in y-direction
            save_path (str): filepath to where the output file should be stored

        Returns: Nothing, it writes a raster file with decreased resolution.

        """
        props = self.dem.GetGeoTransform()
        print('current pixel xsize:', props[1], 'current pixel ysize:', -props[-1])
        options = gdal.WarpOptions(options=['tr'], xRes=new_xres, yRes=new_yres)
        newfile = gdal.Warp(save_path, self.dem, options=options)
        newprops = newfile.GetGeoTransform()
        print('new pixel xsize:', newprops[1], 'new pixel ysize:', -newprops[-1])
        print('file saved in ' + save_path) 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import Warp [as 别名]
def reproject_image(in_raster, out_raster_path, new_projection,  driver = "GTiff",  memory = 2e3, do_post_resample=True):
    """
    Creates a new, reprojected image from in_raster using the gdal.ReprojectImage function.

    Parameters
    ----------
    in_raster
        Either a gdal.Raster object or a path to a raster
    out_raster_path
        The path to the new output raster.
    new_projection
        The new projection in .wkt
    driver
        The format of the output raster.
    memory
        The amount of memory to give to the reprojection.
    do_post_resample
        If set to false, do not resample the image back to the original projection.

    Notes
    -----
    The GDAL reprojection routine changes the size of the pixels by a very small amount; for example, a 10m pixel image
    can become a 10.002m pixel resolution image. To stop alignment issues,
    by deafult this function resamples the images back to their original resolution

    """
    if type(new_projection) is int:
        proj = osr.SpatialReference()
        proj.ImportFromEPSG(new_projection)
        new_projection = proj.ExportToWkt()
    log = logging.getLogger(__name__)
    log.info("Reprojecting {} to {}".format(in_raster, new_projection))
    if type(in_raster) is str:
        in_raster = gdal.Open(in_raster)
    res = in_raster.GetGeoTransform()[1]
    gdal.Warp(out_raster_path, in_raster, dstSRS=new_projection, warpMemoryLimit=memory, format=driver)
    # After warping, image has irregular gt; resample back to previous pixel size
    # TODO: Make this an option
    if do_post_resample:
        resample_image_in_place(out_raster_path, res)
    return out_raster_path 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import Warp [as 别名]
def resample_image_in_place(image_path, new_res):
    """
    Resamples an image in-place using gdalwarp to new_res in metres.
    WARNING: This will make a permanent change to an image! Use with care.

    Parameters
    ----------
    image_path
        Path to the image to be resampled

    new_res
        Pixel edge size in meters

    """
    # I don't like using a second object here, but hey.
    with TemporaryDirectory() as td:
        # Remember this is used for masks, so any averging resample strat will cock things up.
        args = gdal.WarpOptions(
            xRes=new_res,
            yRes=new_res
        )
        temp_image = os.path.join(td, "temp_image.tif")
        gdal.Warp(temp_image, image_path, options=args)

        # Urrrgh. Stupid Windows permissions.
        if sys.platform.startswith("win"):
            os.remove(image_path)
            shutil.copy(temp_image, image_path)
        else:
            shutil.move(temp_image, image_path) 

# 需要导入模块: import gdal [as 别名]
# 或者: from gdal import Warp [as 别名]
def clip_raster(raster_path, aoi_path, out_path, srs_id=4326, flip_x_y = False):
    """
    Clips a raster at raster_path to a shapefile given by aoi_path. Assumes a shapefile only has one polygon.
    Will np.floor() when converting from geo to pixel units and np.absolute() y resolution form geotransform.
    Will also reproject the shapefile to the same projection as the raster if needed.

    Parameters
    ----------
    raster_path
        Path to the raster to be clipped.
    aoi_path
        Path to a shapefile containing a single polygon
    out_path
        Path to a location to save the final output raster

    """
    # https://gis.stackexchange.com/questions/257257/how-to-use-gdal-warp-cutline-option
    with TemporaryDirectory() as td:
        log.info("Clipping {} with {}".format(raster_path, aoi_path))
        raster = gdal.Open(raster_path)
        in_gt = raster.GetGeoTransform()
        srs = osr.SpatialReference()
        srs.ImportFromWkt(raster.GetProjection())
        intersection_path = os.path.join(td, 'intersection')
        aoi = ogr.Open(aoi_path)
        if aoi.GetLayer(0).GetSpatialRef().ExportToWkt() != srs.ExportToWkt():    # Gross string comparison. Might replace with wkb
            log.info("Non-matching projections, reprojecting.")
            aoi = None
            tmp_aoi_path = os.path.join(td, "tmp_aoi.shp")
            reproject_vector(aoi_path, tmp_aoi_path, srs)
            aoi = ogr.Open(tmp_aoi_path)
        intersection = get_aoi_intersection(raster, aoi)
        min_x_geo, max_x_geo, min_y_geo, max_y_geo = intersection.GetEnvelope()
        if flip_x_y:
            min_x_geo, min_y_geo = min_y_geo, min_x_geo
            max_x_geo, max_y_geo = max_y_geo, max_x_geo
        width_pix = int(np.floor(max_x_geo - min_x_geo)/in_gt[1])
        height_pix = int(np.floor(max_y_geo - min_y_geo)/np.absolute(in_gt[5]))
        new_geotransform = (min_x_geo, in_gt[1], 0, max_y_geo, 0, in_gt[5])   # OK, time for hacking
        write_geometry(intersection, intersection_path, srs_id=srs.ExportToWkt())
        clip_spec = gdal.WarpOptions(
            format="GTiff",
            cutlineDSName=intersection_path+r"/geometry.shp",   # TODO: Fix the need for this
            cropToCutline=True,
            width=width_pix,
            height=height_pix,
            dstSRS=srs
        )
        out = gdal.Warp(out_path, raster, options=clip_spec)
        out.SetGeoTransform(new_geotransform)
        out = None 

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