Python osgeo.gdal 模块，ReprojectImage() 实例源码

def reprojectRaster(src,match_ds,dst_filename):
    src_proj = src.GetProjection()
    src_geotrans = src.GetGeoTransform()

    match_proj = match_ds.GetProjection()
    match_geotrans = match_ds.GetGeoTransform()
    wide = match_ds.RasterXSize
    high = match_ds.RasterYSize

    dst = gdal.GetDriverByName('GTiff').Create(dst_filename, wide, high, 1, gdalconst.GDT_Int16)
    dst.SetGeoTransform( match_geotrans )
    dst.SetProjection( match_proj)

    gdal.ReprojectImage(src, dst, src_proj, match_proj, gdalconst.GRA_Bilinear)

    del dst # Flush
    return(gdal.Open(dst_filename,gdalconst.GA_ReadOnly))


#finds the intersecting extent of a series of scenes (left,right,bottom,top are arrays containing the respective lat/lon of the left,right,bottom,top of each image)

def reproject_tiff_from_model(old_name, new_name, model, unit):
    """
    Reprojects an tiff on a tiff model. Can be used to warp tiff.
    Keyword arguments:
    old_name -- the name of the old tiff file
    new_name -- the name of the output tiff file
    model -- the gdal dataset which will be used to warp the tiff
    unit -- the gdal unit in which the operation will be performed
    """
    mem_drv = gdal.GetDriverByName("MEM")

    old = gdal.Open(old_name)

    new = mem_drv.Create(new_name, model.RasterXSize, model.RasterYSize, 1,
                         unit)

    new.SetGeoTransform(model.GetGeoTransform())
    new.SetProjection(model.GetProjection())

    res = gdal.ReprojectImage(old, new, old.GetProjection(),
                              model.GetProjection(), gdal.GRA_NearestNeighbour)

    assert res == 0, 'Error in ReprojectImage'
    arr = new.ReadAsArray()
    new = None
    return arr

def reprojectPanBand(src,match_ds,dst_filename):
    src_proj = src.GetProjection()
    src_geotrans = src.GetGeoTransform()

    match_proj = match_ds.GetProjection()
    match_geotrans = match_ds.GetGeoTransform()
    wide = match_ds.RasterXSize
    high = match_ds.RasterYSize

    dst = gdal.GetDriverByName('GTiff').Create(dst_filename, wide, high, 1, gdalconst.GDT_Int16)
    dst.SetGeoTransform( match_geotrans )
    dst.SetProjection( match_proj)

    gdal.ReprojectImage(src, dst, src_proj, match_proj, gdalconst.GRA_Bilinear)
    return(gdal.Open(dst_filename,gdalconst.GA_ReadOnly))

def reproject_tiff_custom(old_name, new_name, new_x_size, new_y_size,
                          new_geo_transform, new_projection, unit, mode):
    """
    Reprojects an tiff with custom tiff arguments. Can be used to warp tiff.
    If no projection is provided, fallback to old projection.
    Keyword arguments:
    old_name -- the name of the old tiff file
    new_name -- the name of the output tiff file
    new_x_size -- the number of new size in x dimension
    new_y_size -- the number of new size in y dimension
    new_geo_transform -- the new geo transform to apply
    new_projection -- the new projection to use
    unit -- the gdal unit in which the operation will be performed
    mode -- the gdal mode used for warping
    """
    mem_drv = gdal.GetDriverByName("MEM")

    old = gdal.Open(old_name)
    r = old.GetRasterBand(1)
    r.GetNoDataValue()

    # Adds 1 to keep the original zeros (reprojectImage maps NO_DATA to 0)
    old_array = old.ReadAsArray()
    mask = old_array == old.GetRasterBand(1).GetNoDataValue()
    old_array += 1
    old_array[mask] = 0
    temp = mem_drv.Create("temp", old.RasterXSize, old.RasterYSize, 1, unit)
    temp.SetGeoTransform(old.GetGeoTransform())
    temp.SetProjection(old.GetProjection())
    temp.GetRasterBand(1).WriteArray(old_array)

    new = mem_drv.Create(new_name, new_x_size, new_y_size, 1, unit)

    new.SetGeoTransform(new_geo_transform)
    if new_projection is None:
        new.SetProjection(old.GetProjection())
    else:
        new.SetProjection(new_projection)

    res = gdal.ReprojectImage(temp, new, old.GetProjection(), new_projection,
                              mode)

    assert res == 0, 'Error in ReprojectImage'
    arr = new.ReadAsArray()
    mask = arr != 0
    arr -= 1
    arr[~mask] = 0
    new = None
    temp = None

    return arr, mask

def reproject_dataset_example(dataset, dataset_example, method=1):

    # open dataset that must be transformed 
    try:
        if dataset.split('.')[-1] == 'tif':
            g = gdal.Open(dataset)                                  
        else:
            g = dataset    
    except:
            g = dataset            
    epsg_from = Get_epsg(g)    

    # open dataset that is used for transforming the dataset
    try:
        if dataset_example.split('.')[-1] == 'tif':
            gland = gdal.Open(dataset_example)                                  
        else:
            gland = dataset_example      
    except:
            gland = dataset_example              
    epsg_to = Get_epsg(gland)   

    # Set the EPSG codes
    osng = osr.SpatialReference()
    osng.ImportFromEPSG(epsg_to)
    wgs84 = osr.SpatialReference()
    wgs84.ImportFromEPSG(epsg_from)

    # Get shape and geo transform from example              
    geo_land = gland.GetGeoTransform()          
    col=gland.RasterXSize
    rows=gland.RasterYSize

    # Create new raster         
    mem_drv = gdal.GetDriverByName('MEM')
    dest1 = mem_drv.Create('', col, rows, 1, gdal.GDT_Float32)
    dest1.SetGeoTransform(geo_land)
    dest1.SetProjection(osng.ExportToWkt())

    # Perform the projection/resampling
    if method is 1:             
        gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_NearestNeighbour)
    if method is 2:             
        gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_Bilinear)
    if method is 3:             
        gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_Lanczos)
    if method is 4:             
        gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_Average)
    return(dest1)