Vorrei creare alcuni raster di elevazione/altezza utilizzando python, gdal e numpy. Sono bloccato sul NumPy (. E, probabilmente, pitone e GDAL)creazione di altezza/altezza campo gdal numpy python
In NumPy, sto tentando il seguente:
>>> 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
da OSGeo importazione gdal da 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 figura mi manca qualcosa di semplice e attendo con ansia il vostro consiglio.
Grazie,
Chris
(continua dopo)
- terragendataset.cpp, v 1,2 *
- Progetto: Terragen (tm) TER driver
- Scopo: Reader per documenti Terragen TER
- Autore: Ray Gardener, Daylon Graphics Ltd. *
- Parti di questo modulo derivato dal driver GDAL da
- Frank Warmerdam, vedere http://www.gdal.org
mie scuse in anticipo per Ray e Frank Gardener Warmerdam.
Terragen note formato:
Per la scrittura: SCAL = distanza in metri gridpost hv_px = hv_m/SCAL span_px = span_m/SCAL offset = vedere TerragenDataset :: write_header() scala = vedi TerragenDataset: : write_header() hv fisico = (hv_px - offset) * 65.536,0/scala
Diciamo ai chiamanti che:
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.
Questo mi dice che prima del mio sopra WriteArray (somearray) devo impostare sia la trasformazione geografica e SetProjection e SetUnitType cosa c'è da lavorare (potenzialmente)
dal tutorial GDAL API: Python importazione OSR importazione 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
Le mie scuse per la creazione di un troppo lungo post e una confessione. Se e quando avrò ragione, sarebbe bello avere tutte le note in un posto (il post lungo).
The Confession:
precedenza ho scattato una foto (jpeg), convertito in un GeoTIFF, ed importato come piastrelle in un database PostGIS. Sto ora tentando di creare un raster di elevazione su cui applicare l'immagine. Questo probabilmente sembra sciocco, ma c'è un artista che desidero onorare, mentre allo stesso tempo non offendo coloro che hanno lavorato assiduamente per creare e migliorare questi grandi strumenti.
L'artista è belga, quindi i metri sarebbero in ordine. Lavora nella Lower Manhattan, NY quindi, UTM 18. Ora alcuni SetGeoTransform sono ragionevoli. L'immagine sopra menzionata è w = 3649/h = 2736, dovrò riflettere su questo.
un altro tentativo:
>>> 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
Chiaramente sempre più vicino, ma chiaro se il SetUTM (18,1) è stato preso. Si tratta di un 4x4 nel fiume Hudson o in un Local_CS (sistema di coordinate)? Cos'è un fallito silenzioso?
lettura più
// 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 metri è un grazioso piccolo arco di logica.
Quindi, forse questo è buono come si arriva. SetGeoTransform ottiene le unità corrette, imposta la scala e hai il tuo Terragen World Space.
Pensiero finale: non programmo, ma in parte posso seguire. Detto questo, ho un pò chiesto prima se e poi che i dati sembravano nel mio piccolo Terragen World Space (le seguenti grazie alla http://www.gis.usu.edu/~chrisg/python/2009/ settimana 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)
>>>
Quindi questo è gratificante. Immagino la differenza tra il numpy usato sopra c e questo risultato va alle azioni di applicare Float16 attraverso questo molto piccolo intervallo logico .
E a '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
>>>
quasi completamente gratificante, anche se sembra come se fossi in La Concordia Perù. Quindi ho per capire come dire - come più a nord, come New York North. SetUTM prende un terzo elemento che suggerisce "quanto lontano" nord o sud. Sembra che ieri ho trovato un grafico UTM con zone di lettere, qualcosa come C all'equatore e forse T o S per l'area di New York.
In realtà ho pensato che SetGeoTransform stava essenzialmente stabilendo il nord in alto a sinistra e ad est e questo stava influenzando la parte di "quanto lontano nord/sud" di SetUTM. Vai a gdal-dev.
E poi ancora:
Paddington Orso è andato in Perù perché aveva un biglietto. Sono arrivato perché ho detto che era dove volevo essere. Terragen, lavorando come fa, mi ha dato il mio spazio pixel. Le successive chiamate a srs hanno agito su hf2 (SetUTM), ma est e nord sono stati stabiliti con Terragen, quindi l'UTM 18 è stato impostato, ma in un riquadro di delimitazione all'equatore. Abbastanza buono.
gdal_translate mi ha portato a New York. Sono in Windows quindi una riga di comando.e il risultato;
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>
Quindi, siamo di nuovo a New York. Ci sono probabilmente modi migliori per affrontare tutto questo. Io dovevo avere un obiettivo che accettava Crea mentre stavo postando/improvvisando un set di dati da Numpy. Ho bisogno di guardare altri formati che consentono di creare. L'elevazione in GeoTiff è una possibilità (credo)
I miei ringraziamenti a tutti i commenti, suggerimenti e gentili suggerimenti per una lettura appropriata. Creare mappe in Python è divertente!
Chris
Qual è la domanda? Stai scrivendo degli zeri su un file ... Che cosa vuoi fare? La scrittura degli zeri non funziona? Se si desidera scrivere la matrice numpy nel file, passarla invece della matrice di zeri che si sta creando. (Potrebbe essere necessario passare in 'c.astype (numpy.float32)', se si sta creando una banda a 32 bit e 'c' è un array a 64 bit (che dipenderà dalla piattaforma in cui ci si trova)). –
[Il driver] (http://www.gdal.org/frmt_terragen.html) supporta solo 'gdal.GDT_Int16'-non float32 –
@JoeKingston - gli zeri funzionavano, ma volevo passare c come float 32, come Sto creando i miei dati di altezza. – Chris