CODEWITHSUNDEEP
pythonscipygisgdal
jreid28
jreid28

Reputation: 21

GDAL Affine Coefficients from Lat/Lon Meshgrids in Python

I am having some problems with Affine tranformation coefficients while creating a new GeoTIFF file. What I am doing is ETL on a scientific dataset that results in a 2D Ndarray along with a set of meshgrid Ndarrays that contain Lat and Lon. Both the meshgrids and the dataset arrays have the same dimensions of 645 x 980. From what I understand the GeoTIFF requires a list of Affine coefficients when created from Python GDAL via the SetGeoTransform() method. The list has the form of [xllcorner, xrotation, x_cellsize, yllcorner, yrotation, y_cellsize]. My approach to this is similar to what is outlined here: http://adventuresindevelopment.blogspot.com/2008/12/python-gdal-adding-geotiff-meta-data.html

At this point is where I run into problems. I calculate the xllcorner and the yllcorner using the min() method for the two meshgrid arrays for lat & lon respectively, and I manually calculate the x and y cellsize by applying the formula [max-min]/dimension size with the x dimension size being the x axis size for the lons meshgrid and the y dimension size being the y axis size for the lats meshgrid. When I apply this method and try to write out the array band via GetRasterBand().WriteArray() I get this error message:

Traceback (most recent call last):
    ...
    raise ValueError("array larger than output file, or offset off edge")
ValueError: array larger than output file, or offset off edge

Therefore I assume that I have composed my affine coefficients incorrectly but given the data this makes no sense to me. I even made sure that the Spatial Reference System was set to WGS:84 before attempting the affine coefficient creation. So my question is how to properly create the Affine coefficients with lat/lon meshgrids and a data array that share common dimensions? I think my cell size calculation can't simply be lat/lon differences; but I am not sure.

Upvotes: 2

Views: 1032

Answers (1)

Mike T
Mike T

Reputation: 43682

This error is typically shown when the expected array shape does not match. For instance, see what shape the expected shape is with:

band = src.GetRasterBand(1)
arr = band.ReadAsArray()
print(arr.shape)  # (656L, 515L)

This will need to be the shape of the numpy array to be written:

assert other_array.shape == arr.shape
band.WriteArray(other_array)

And to raise the same ValueError, change the shape so it is longer in one dimension, e.g.:

band.WriteArray(other_array.T)

As for affine transformations, this is probably not raising any errors, as it is often just stored as data. GIS rasters typically register the world coordinate in the upper-left corner, and use a -dy value to count rows downwards. However, using a lower-left corner with +dy is usually fine by most software. It will just be upside down when comparing the array as a printed matrix versus mapped raster.

Upvotes: 2

Related Questions