Python GDAL:使用其他文件进行投影的地理参考数组 [英] Python GDAL: Georeference array using other file for projection

问题描述

我有一个数据数组，对于每个点，我都知道该点的纬度和经度，因此我想将数据从另一个文件中提取出来，然后将其写入GTiff.如何正确地理参考新文件?

I have an array of data, for each point I know the latitude and longitude of that point, and I'd like to write the data to a GTiff with projection taken from another file. How do I properly georeference the new file?

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

This is what I'm attempting just now:

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格式)您的数据.除非您要进行某种形式的重新投影或插值，否则不需要目标栅格".

If all you want to do is save the data to a raster for use in QGIS, you can simply construct a new Geotiff (or any other GDAL format) from your data. There is no need for a 'target raster' unless you want to do some form of reprojection or interpolation.

这里是一个例子:

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处理边缘，因此有必要增加一半像素大小.

In this example i assumed that your lat/lon's refer to the center of a pixel, since GDAL works with the edge, adding half a pixelsize is necessary.

这篇关于Python GDAL:使用其他文件进行投影的地理参考数组的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！