raster_conversions.py 文件源码

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)