Wed 05 September 2007
I've been working alot in Surfer these days; an excellent geostats and surface mapping package. I was very happy to find that GDAL read it's .grd binary format until I noticed the output from gdalinfo:
> C:\Workspace\Temp\interpolation>gdalinfo svpce_5.grd
Driver: GS7BG/Golden Software 7 Binary Grid (.grd)
Files: svpce_5.grd
Size is 555, 339
Coordinate System is `'
Origin = (383371.000000000000000,3764907.000000000000000)
Pixel Size = (0.500000000000000,0.500000000000000)
Corner Coordinates:
Upper Left ( 383371.000, 3764907.000)
Lower Left ( 383371.000, 3765076.500)
Upper Right ( 383648.500, 3764907.000)
Lower Right ( 383648.500, 3765076.500)
Center ( 383509.750, 3764991.750)
Band 1 Block=555x1 Type=Float64, ColorInterp=Undefined
NoData Value=1.70141e+038
Notice that upper Y value is south of the lower Y value! Basically the raster lines order is reversed (bottom-to-top instead of the normal raster orientation of top-to-bottom). I've also experienced the same issue with some NetCDF files so I thought it would be good to have a generic solution to the problem.
So I hacked up the gdal_merge.py script (distributed with gdal, fwtools, etc) and created a raster flip script that will invert the image along the y axis and retain the georeferencing and metadata. The resulting flip_raster.py script seems to work pretty well though it is far from tested.
Here's an example:
The standard gdal_translate method (which doesn't account for the inverted coordinate space):
gdal_translate -of GTiff krig1.grd krig1_translate.tif
And the flipped raster method:
flip_raster.py -o krig1_flip.tif -of GTiff krig1.grd
And we're good. gdalinfo confirms that we have the same extents, pixel sizes, metadata, etc as the original dataset.