Tengo una matriz de datos, y para cada punto de datos sé la latitud y la longitud. Me gustaría guardarlo como un GTiff con la misma proyección que otros rásters que tengo. Esto es lo que he intentado hasta ahora, pero no tuve suerte.
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)
Pero esto da como resultado el siguiente mensaje de error:
File "TES2GtifBBB.py", line 25, in CreateGeoTiff
DataSet.GetRasterBand(1).WriteArray( Array )
File "/Library/Frameworks/GDAL.framework/Versions/1.9/Python/2.7/site packages/osgeo/gdal.py", line 1082, in WriteArray
return gdalnumeric.BandWriteArray( self, array, xoff, yoff )
File "/Library/Frameworks/GDAL.framework/Versions/1.9/Python/2.7/site packages/osgeo/gdal_array.py", line 256, in BandWriteArray
raise ValueError("array larger than output file, or offset off edge")
ValueError: array larger than output file, or offset off edge