Python GDAL：使用其他文件进行投影的地理参考数组答案

【问题标题】：Python GDAL: Georeference array using other file for projectionPython GDAL：使用其他文件进行投影的地理参考数组
【发布时间】：2013-04-19 01:15:07
【问题描述】：

我有一个数据数组，对于每个点，我都知道该点的纬度和经度，并且我想将数据写入 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)

【问题讨论】：

你想做什么？如果您使用SourceProjection = TargetProjection.CloneGeogCS()，您是否只想将数据值插入到目标网格中？使用 osr.CoordinateTransformation 重新投影对于栅格/网格没有用处，因为您的源网格和目标网格可能没有“对齐”，因此如果您想要填充目标网格，则需要进行某种形式的插值。
@RutgerKassies 是的，我知道重新投影可能会导致某种形式的插值，但在一个可能不会太关键的小区域内。我只想将我的数组保存在一个表单中，以便我可以在 QGIS 中将它与我的其他栅格文件一起打开。
我把它写成另一个问题 (gis.stackexchange.com/questions/58517/…)，然后由于某种原因这个问题从我的个人资料中消失了，所以我认为我没有点击保存或其他东西并在这里重写它。现在原件出现了，有一个答案，也许有必要重复的人可以将此问题标记为重复？

标签： python gdal

【解决方案1】：

如果您只想将数据保存到栅格以在 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 与边缘一起工作，所以需要添加半个像素大小。

【讨论】：