creating elevation/height field gdal numpy python

Question

I would like to create some elevation/heightfield rasters using python, gdal and numpy. I'm stuck on numpy (and probably python and gdal.)

In numpy, I've been attempting the following:

>>> a= numpy.linspace(4,1,4, endpoint=True)
>>> b= numpy.vstack(a)
>>> c=numpy.repeat(b,4,axis=1)
>>> c
array([[ 4.,  4.,  4.,  4.],
       [ 3.,  3.,  3.,  3.],
       [ 2.,  2.,  2.,  2.],
       [ 1.,  1.,  1.,  1.]]) #This is the elevation data I want

from osgeo import gdal from gdalconst import *

>>> format = "Terragen"
>>> driver = gdal.GetDriverByName(format)    
>>> dst_ds = driver.Create('test.ter', 4,4,1,gdal.GDT_Float32, ['MINUSERPIXELVALUE = 1',                                                             'MAXUSERPIXELVALUE= 4']) 
>>> raster = numpy.zeros((4,4), dtype=numpy.float32) #this is where I'm messing up
>>> dst_ds.GetRasterBand(1).WriteArray(raster) # gives the null elevation data I asked for in (raster) 
0
>>> dst_ds = None

I figure I'm missing something simple and look forward to your advice.

Thanks,

Chris

(continued later)

• terragendataset.cpp,v 1.2 *
  • Project: Terragen(tm) TER Driver
  • Purpose: Reader for Terragen TER documents
  • Author: Ray Gardener, Daylon Graphics Ltd. *
  • Portions of this module derived from GDAL drivers by
  • Frank Warmerdam, see http://www.gdal.org

My apologies in advance to Ray Gardener and Frank Warmerdam.

Terragen format notes:

For writing: SCAL = gridpost distance in meters hv_px = hv_m / SCAL span_px = span_m / SCAL offset = see TerragenDataset::write_header() scale = see TerragenDataset::write_header() physical hv = (hv_px - offset) * 65536.0/scale

We tell callers that:

    Elevations are Int16 when reading, 
    and Float32 when writing. We need logical 
    elevations when writing so that we can 
    encode them with as much precision as possible 
    when going down to physical 16-bit ints.
    Implementing band::SetScale/SetOffset won't work because 
    it requires callers to know format write details.
    So we've added two Create() options that let the 
    caller tell us the span's logical extent, and with 
    those two values we can convert to physical pixels.

            ds::SetGeoTransform() lets us establish the
        size of ground pixels.
    ds::SetProjection() lets us establish what
        units ground measures are in (also needed 
        to calc the size of ground pixels).
    band::SetUnitType() tells us what units
        the given Float32 elevations are in.

This tells me that prior to my above WriteArray(somearray) I have to set both the GeoTransform and SetProjection and SetUnitType for things to work (potentially)

From the GDAL API Tutorial: Python import osr import numpy

dst_ds.SetGeoTransform( [ 444720, 30, 0, 3751320, 0, -30 ] )

srs = osr.SpatialReference()
srs.SetUTM( 11, 1 )
srs.SetWellKnownGeogCS( 'NAD27' )
dst_ds.SetProjection( srs.ExportToWkt() )

raster = numpy.zeros( (512, 512), dtype=numpy.uint8 ) #we've seen this before   
dst_ds.GetRasterBand(1).WriteArray( raster )

# Once we're done, close properly the dataset
dst_ds = None 

My apologies for creating an overly long post and a confession. If and when I get this right, it would be nice to have all the notes in one place (the long post).

The Confession:

I have previously taken a picture (jpeg), converted it to a geotiff, and imported it as tiles into a PostGIS database. I am now attempting to create elevation raster over which to drape the picture. This probably seems silly, but there is an artist whom I wish to honor, while at the same time not offending those who have worked assiduously to create and improve these great tools.

The artist is Belgian so meters would be in order. She works in lower Manhattan, NY so, UTM 18. Now some reasonable SetGeoTransform. The above mentioned picture is w=3649/h=2736, I'll have to give some thought to this.

Another try:

>>> format = "Terragen"
>>> driver = gdal.GetDriverByName(format)
>>> dst_ds = driver.Create('test_3.ter', 4,4,1, gdal.GDT_Float32, ['MINUSERPIXELVALUE=1','MAXUSERPIXELVALUE-4']) 
>>> type(dst_ds)
<class 'osgeo.gdal.Dataset'>
>>> import osr
>>> dst_ds.SetGeoTransform([582495, 1, 0.5, 4512717, 0.5, -1]) #x-res 0.5, y_res 0.5 a guess
0
>>> type(dst_ds)
<class 'osgeo.gdal.Dataset'>
>>> srs = osr.SpatialReference()
>>> srs.SetUTM(18,1)
0
>>> srs.SetWellKnownGeogCS('WGS84')
0
>>> dst_ds.SetProjection(srs.ExportToWkt())
0
>>> raster = c.astype(numpy.float32)
>>> dst_ds.GetRasterBand(1).WriteArray(raster)
0
>>> dst_ds = None
>>> from osgeo import gdalinfo
>>> gdalinfo.main(['foo', 'test_3.ter'])
Driver: Terragen/Terragen heightfield
Files: test_3.ter
Size is 4, 4
Coordinate System is:
LOCAL_CS["Terragen world space",
    UNIT["m",1]]
Origin = (0.000000000000000,0.000000000000000)
Pixel Size = (1.000000000000000,1.000000000000000)
Metadata:
  AREA_OR_POINT=Point
Corner Coordinates:
Upper Left  (   0.0000000,   0.0000000) 
Lower Left  (   0.0000000,   4.0000000) 
Upper Right (   4.0000000,   0.0000000) 
Lower Right (   4.0000000,   4.0000000) 
Center      (   2.0000000,   2.0000000) 
Band 1 Block=4x1 Type=Int16, ColorInterp=Undefined
  Unit Type: m
Offset: 2,   Scale:7.62939453125e-05
0

Clearly getting closer but unclear if the SetUTM(18,1) was picked up. Is this a 4x4 in the Hudson River or a Local_CS(coordinate system)? What's a silent fail?

More Reading

// Terragen files aren't really georeferenced, but 
// we should get the projection's linear units so 
// that we can scale elevations correctly.

// Increase the heightscale until the physical span
// fits within a 16-bit range. The smaller the logical span,
// the more necessary this becomes.

4x4 meters is a pretty small logical span.

So, perhaps this is as good as it gets. The SetGeoTransform gets the units right, sets the scale and you have your Terragen World Space.

Final Thought: I don't program, but to an extent I can follow along. That said, I kinda wondered first if and then what the data looked like in my little Terragen World Space (the following thanks to http://www.gis.usu.edu/~chrisg/python/2009/ week 4):

>>> fn = 'test_3.ter'
>>> ds = gdal.Open(fn, GA_ReadOnly)
>>> band = ds.GetRasterBand(1)
>>> data = band.ReadAsArray(0,0,1,1)
>>> data
array([[26214]], dtype=int16)
>>> data = band.ReadAsArray(0,0,4,4)
>>> data
array([[ 26214,  26214,  26214,  26214],
       [ 13107,  13107,  13107,  13107],
       [     0,      0,      0,      0],
       [-13107, -13107, -13107, -13107]], dtype=int16)
>>>

So this is gratifying. I imagine the difference between the above used numpy c and this result goes to the actions of applying Float16 across this very small logical span.

And on to 'hf2':

>>> src_ds = gdal.Open('test_3.ter')
>>> dst_ds = driver.CreateCopy('test_6.hf2', src_ds, 0)
>>> dst_ds.SetGeoTransform([582495,1,0.5,4512717,0.5,-1])
0
>>> srs = osr.SpatialReference()
>>> srs.SetUTM(18,1)
0
>>> srs.SetWellKnownGeogCS('WGS84')
0
>>> dst_ds.SetProjection( srs.ExportToWkt())
0
>>> dst_ds = None
>>> src_ds = None
>>> gdalinfo.main(['foo','test_6.hf2'])
Driver: HF2/HF2/HFZ heightfield raster
Files: test_6.hf2
   test_6.hf2.aux.xml
Size is 4, 4
Coordinate System is:
PROJCS["UTM Zone 18, Northern Hemisphere",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        TOWGS84[0,0,0,0,0,0,0],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9108"]],
    AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-75],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["Meter",1]]
Origin = (0.000000000000000,0.000000000000000)
Pixel Size = (1.000000000000000,1.000000000000000)
Metadata:
 VERTICAL_PRECISION=1.000000
Corner Coordinates:
Upper Left  (   0.0000000,   0.0000000) ( 79d29'19.48"W,  0d 0' 0.01"N)
Lower Left  (   0.0000000,   4.0000000) ( 79d29'19.48"W,  0d 0' 0.13"N)
Upper Right (   4.0000000,   0.0000000) ( 79d29'19.35"W,  0d 0' 0.01"N)
Lower Right (   4.0000000,   4.0000000) ( 79d29'19.35"W,  0d 0' 0.13"N)
Center      (   2.0000000,   2.0000000) ( 79d29'19.41"W,  0d 0' 0.06"N)
Band 1 Block=256x1 Type=Float32, ColorInterp=Undefined
Unit Type: m
0
>>> 

Nearly completely gratifying, though it looks like I'm in La Concordia Peru. So I have to figure out how to say- like more north, like New York North. Does SetUTM take a third element suggesting 'how far' north or south. Seems I ran across a UTM chart yesterday that had letter label zones, something like C at the equator and maybe T or S for New York area.

I actually thought the SetGeoTransform was essentially establishing the top left northing and easting and this was influencing the that 'how far north/south' part of SetUTM. Off to gdal-dev.

And later still:

Paddington Bear went to Peru because he had a ticket. I got there because I said that's where I wanted to be. Terragen, working as it does, gave me my pixel space. The subsequent calls to srs acted upon the hf2 (SetUTM), but the easting and northing were established under Terragen so the UTM 18 got set, but in a bounding box at the equator. Good enough.

gdal_translate got me to New York. I'm in windows so a command line. and the result;

C:\Program Files\GDAL>gdalinfo c:/python27/test_9.hf2
Driver: HF2/HF2/HFZ heightfield raster
Files: c:/python27/test_9.hf2
Size is 4, 4
Coordinate System is:
PROJCS["UTM Zone 18, Northern Hemisphere",
GEOGCS["NAD83",
    DATUM["North_American_Datum_1983",
        SPHEROID["GRS 1980",6378137,298.257222101,
            AUTHORITY["EPSG","7019"]],
        TOWGS84[0,0,0,0,0,0,0],
        AUTHORITY["EPSG","6269"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4269"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-75],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["Meter",1]]
Origin = (583862.000000000000000,4510893.000000000000000)
Pixel Size = (-1.000000000000000,-1.000000000000000)
Metadata:
VERTICAL_PRECISION=0.010000
Corner Coordinates:
Upper Left  (  583862.000, 4510893.000) ( 74d 0'24.04"W, 40d44'40.97"N)
Lower Left  (  583862.000, 4510889.000) ( 74d 0'24.04"W, 40d44'40.84"N)
Upper Right (  583858.000, 4510893.000) ( 74d 0'24.21"W, 40d44'40.97"N)
Lower Right (  583858.000, 4510889.000) ( 74d 0'24.21"W, 40d44'40.84"N)
Center      (  583860.000, 4510891.000) ( 74d 0'24.13"W, 40d44'40.91"N)
Band 1 Block=256x1 Type=Float32, ColorInterp=Undefined
Unit Type: m

C:\Program Files\GDAL>

So, we're back in NY. There are probably better ways to approach all this. I had to have a target that accepted Create as I was postulating/improvising a dataset from numpy as well. I need to look at other formats that allow create. Elevation in GeoTiff is a possibility (I think.)

My thanks to all Comments, suggestions, and gentle nudges toward appropriate reading. Making maps in python is fun!

Chris

Answer 1

You aren't too far off. Your only problem is the syntax issues with the GDAL creation options. Replace:

['MINUSERPIXELVALUE = 1','MAXUSERPIXELVALUE= 4']

with no spaces before or after the key/value pairs:

['MINUSERPIXELVALUE=1', 'MAXUSERPIXELVALUE=4']

Check type(dst_ds) and it should be <class 'osgeo.gdal.Dataset'> rather than <type 'NoneType'> for a silent fail as was the case for the above.

---

Update By default, warnings or exceptions are not shown. You can enable them via gdal.UseExceptions() to see what's ticking underneath, e.g.:

>>> from osgeo import gdal
>>> gdal.UseExceptions()
>>> driver = gdal.GetDriverByName('Terragen')
>>> dst_ds = driver.Create('test.ter', 4,4,1,gdal.GDT_Float32, ['MINUSERPIXELVALUE = 1', 'MAXUSERPIXELVALUE= 4'])
Warning 6: Driver Terragen does not support MINUSERPIXELVALUE  creation option
>>> dst_ds = driver.Create('test.ter', 4,4,1,gdal.GDT_Float32, ['MINUSERPIXELVALUE=1', 'MAXUSERPIXELVALUE=4'])
>>> type(dst_ds)
<class 'osgeo.gdal.Dataset'>