Python GDAL : Georeference array using other file for projection

我有一个数据数组，对于每个点，我都知道该点的纬度和经度，我想将数据写入 GTiff，并使用从另一个文件获取的投影。如何正确地对新文件进行地理配准？

这就是我现在正在尝试的:

import numpy as np 
import gdal 
from gdalconst import * 
from osgeo import osr 
 
def GetGeoInfo(FileName): 
    SourceDS = gdal.Open(FileName, GA_ReadOnly) 
    GeoT = SourceDS.GetGeoTransform() 
    Projection = osr.SpatialReference() 
    Projection.ImportFromWkt(SourceDS.GetProjectionRef()) 
    return GeoT, Projection 
 
def CreateGeoTiff(Name, Array, driver,  
                  xsize, ysize, GeoT, Projection): 
    DataType = gdal.GDT_Float32 
    NewFileName = Name+'.tif' 
    # Set up the dataset 
    DataSet = driver.Create( NewFileName, xsize, ysize, 1, DataType ) 
            # the '1' is for band 1. 
    DataSet.SetGeoTransform(GeoT) 
    DataSet.SetProjection( Projection.ExportToWkt() ) 
    # Write the array 
    DataSet.GetRasterBand(1).WriteArray( Array ) 
    return NewFileName 
 
def ReprojectCoords(x, y,src_srs,tgt_srs): 
    trans_coords=[] 
    transform = osr.CoordinateTransformation( src_srs, tgt_srs) 
    x,y,z = transform.TransformPoint(x, y) 
    return x, y 
 
# Some Data 
Data = np.random.rand(5,6) 
Lats = np.array([-5.5, -5.0, -4.5, -4.0, -3.5]) 
Lons = np.array([135.0, 135.5, 136.0, 136.5, 137.0, 137.5]) 
 
# A raster file that exists in the same approximate aregion. 
RASTER_FN = 'some_raster.tif' 
 
# Open the raster file and get the projection, that's the 
# projection I'd like my new raster to have, it's 'projected', 
# i.e. x, y values are numbers of pixels. 
GeoT, TargetProjection, DataType = GetGeoInfo(RASTER_FN) 
# Meanwhile my raster is currently in geographic coordinates. 
SourceProjection = TargetProjection.CloneGeogCS() 
 
# Get the corner coordinates of my array 
LatSize, LonSize = len(Lats), len(Lons) 
LatLow, LatHigh = Lats[0], Lats[-1] 
LonLow, LonHigh = Lons[0], Lons[-1] 
# Reproject the corner coordinates from geographic 
# to projected... 
TopLeft = ReprojectCoords(LonLow, LatHigh, SourceProjection, TargetProjection) 
BottomLeft = ReprojectCoords(LonLow, LatLow, SourceProjection, TargetProjection) 
TopRight = ReprojectCoords(LonHigh, LatHigh, SourceProjection, TargetProjection) 
# And define my Geotransform 
GeoTNew = [TopLeft[0],  (TopLeft[0]-TopRight[0])/(LonSize-1), 0, 
           TopLeft[1], 0, (TopLeft[1]-BottomLeft[1])/(LatSize-1)] 
 
# I want a GTiff 
driver = gdal.GetDriverByName('GTiff') 
# Create the new file... 
NewFileName = CreateGeoTiff('Output', Data, driver, LatSize, LonSize, GeoTNew, TargetProjection) 

请您参考如下方法:

如果您只想将数据保存到栅格中以便在 QGIS 中使用，您可以简单地从您的数据构建一个新的 Geotiff(或任何其他 GDAL 格式)。除非您想进行某种形式的重投影或插值，否则不需要“目标栅格”。

这是一个例子:

import gdal 
import osr 
import numpy as np 
 
data = np.random.rand(5,6) 
lats = np.array([-5.5, -5.0, -4.5, -4.0, -3.5]) 
lons = np.array([135.0, 135.5, 136.0, 136.5, 137.0, 137.5]) 
 
xres = lons[1] - lons[0] 
yres = lats[1] - lats[0] 
 
ysize = len(lats) 
xsize = len(lons) 
 
ulx = lons[0] - (xres / 2.) 
uly = lats[-1] - (yres / 2.) 
 
driver = gdal.GetDriverByName('GTiff') 
ds = driver.Create('D:\\test.tif', xsize, ysize, 1, gdal.GDT_Float32) 
 
# this assumes the projection is Geographic lat/lon WGS 84 
srs = osr.SpatialReference() 
srs.ImportFromEPSG(4326) 
ds.SetProjection(srs.ExportToWkt()) 
 
gt = [ulx, xres, 0, uly, 0, yres ] 
ds.SetGeoTransform(gt) 
 
outband = ds.GetRasterBand(1) 
outband.WriteArray(data) 
 
ds = None 

在这个例子中，我假设你的纬度/经度是指像素的中心，因为 GDAL 使用边缘，所以需要添加一半的像素大小。