I would like to create some raster / elevation fields using python, gdal and numpy. I am stuck on numpy (and probably python and gdal.)
In numpy, I am trying to do 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.]])
from osgeo import gdal from import gdalconst *
>>> 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)
I suppose I missed something simple and look forward to your advice.
Thanks,
Chris
(continued later)
- terragendataset.cpp, v 1.2 *
- Project: TER Terragen (tm) driver
- Purpose: reading TER Terragen documents
- Posted by Ray Gardener, Daylon Graphics Ltd. *
- Parts of this module derived from GDAL drivers
- Frank Warmerdam, see http://www.gdal.org
My apologies in advance to Ray Gardner and Frank Warmerdam.
Terragen format notes:
To write: SCAL = grid 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 subscribers 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 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 before mine above WriteArray (somearray) I need to set both GeoTransform and SetProjection and SetUnitType 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 )
I apologize for creating too long a post and confession. If and when I get this right, it would be nice to have all the notes in one place (long post).
Confession:
Earlier, I took a snapshot (jpeg), converted it to a geotype and imported it as tiles into a PostGIS database. Now I'm trying to create a bitmap bitmap on which to drape the image. It may seem silly, but there is an artist whom I want to honor, but at the same time not insulting those who have worked hard to create and improve these wonderful tools.
The artist is Belgian, so the counters will be fine. She works in lower Manhattan, NY, so UTM 18. Now some reasonable SetGeoTransform. The above image is w = 3649 / h = 2736, I will have to think about it.
Another attempt:
>>> 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])
Obviously, an approximation, but it is not clear whether SetUTM (18.1) has been established. Is it 4x4 in the Hudson River or Local_CS (coordinate system)? What is silent failure?
Read More
4x4 meters is a fairly small logical range.
So perhaps this is as good as it gets. SetGeoTransform gets the units correctly, sets the scale, and you have your terragen terrestrial space.
Final thought: I am not programming, but to the extent that I can follow. However, at first I wondered if and what then the data looked like in my small outer space Terragen (below 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 it pleases. I believe that the difference between the aforementioned use of numpy c and this result relates to the actions of applying Float16 over this very small logical interval.
And on "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 >>>
Almost completely satisfied, although it seems like I'm in La Concordia Peru. So I have to understand how to say how much more in the north, for example, New York North. Does SetUTM include a third element suggesting "how far" north or south. I seem to have come across a UTM diagram yesterday that had zones with letter marks, something like C at the equator and possibly T or S for the New York area.
I really thought that SetGeoTransform essentially sets the upper left north and east, and this affected how far north / south goes into SetUTM. Off Gdal-dev.
And later still:
The Paddington Bear went to Peru because he had a ticket. I got there because I said where I want. Terragen, working like him, gave me his space in pixels. Subsequent srs calls acted on hf2 (SetUTM), but east and north were set under Terragen, so UTM 18 was installed, but in a bounding box at the equator. Good enough.
gdal_translate brought me to New York. I am in windows so command line. and 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 are back in New York. There are probably better ways to get closer to all of this. I had to have a goal that accepted Create, as I postulated / improvised a data set from numpy. I need to look at other formats that allow me to create. Elevation in GeoTiff is an opportunity (I think.)
My thanks to all the comments, suggestions, and gentle pushings for appropriate reading. Creating maps in python is fun!
Chris