Python GDAL : 다른 파일을 투영에 사용하는 지형 참조 배열

Question

각 포인트에 대해 해당 포인트의 위도와 경도를 알고 있고 다른 데이터에서 투영 된 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) 

Answer 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가 반 pixelsize가 필요한 추가 에지에서 작동 이후의은, 화소의 중심을 참조.