的Python GDAL：地理配準陣列使用其它文件用於投影

我有數據的陣列，對於每個點我知道點的緯度和經度，以及我想將數據寫入到與投影一個GTiff從另一個取文件。我如何正確地對新文件進行地理配準？的Python GDAL：地理配準陣列使用其它文件用於投影

這就是我試圖剛纔：

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)

來源

2013-04-19 EddyTheB

你到底想幹什麼？如果你使用'SourceProjection = TargetProjection.CloneGeogCS（）'，不要你只是想你的數據值插值到您的目標格？用osr.CoordinateTransformation心不是有用的柵格重新投影/柵板，因爲你的源和目標網格也有可能不是「對齊」的，所以需要，如果你想填充目標網格做一些形式的內插。 –

@RutgerKassies是的，我知道，重新投影可能會導致某種形式的內插，但在小面積上，可能不會太挑剔。我只想將我的數組保存在一個表單中，以便我可以使用其他光柵文件在QGIS中打開它。 – EddyTheB

我寫這是另一個問題（http://gis.stackexchange.com/questions/58517/python-gdal-save-array-as-raster-with-projection-from-other-file/58551?noredirect=1# 58551），然後出於某種原因，這個問題從我的個人資料中消失了，所以我認爲我沒有點擊保存或其他東西，並在這裏重寫。現在原來的答案已經出來了，有一個答案，或許某個具有必要條件的人可以將這個問題標記爲重複的東西？ – EddyTheB

如果你想要做的就是將數據保存到在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的工作原理與優勢，加入半個pixelsize是必要的。

來源

2013-04-22 11:28:46

的Python GDAL：地理配準陣列使用其它文件用於投影

回答

相關問題