2

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

1 回答 1

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 与边缘一起工作，因此需要添加半个像素大小。

于 2013-04-22T11:28:46.493 回答