tworzenie pola wysokości/wysokości gdal numpy python

Question

Chciałbym utworzyć rastry elewacji/wysokości za pomocą python, gdal i numpy. Utknąłem na numpy (. I prawdopodobnie Python i gdal)

W numpy, byłem próbując następujące:

>>> 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 

z OSGeo importowej gdal z 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 

Myślę, że brakuje mi czegoś prostego i czekam na twoją radę.

Dzięki,

Chris

(kontynuowane później)

• terragendataset.cpp, v 1,2 *
  • Projekt: Terragen (tm) TER Kierowca
  • Cel: Czytnik dokumentów Terragen TER
  • Autor: Ray Gardener, Daylon Graphics Ltd. *
  • Fragmenty tego modułu pochodzącego z kierowców gdal przez
  • Frank Warmerdam patrz http://www.gdal.org

moje przeprosiny z wyprzedzeniem Ray Ogrodnik i Frank Warmerdam.

Terragen formacie Uwagi:

do zapisu: SCAL = odległość w metrach gridpost hv_px = hv_m/SCAL span_px = span_m/SCAL offset = patrz TerragenDataset :: write_header() skala = patrz TerragenDataset: : write_header() fizyczny hv = (hv_px - offset) * 65536,0/skala

Powiemy rozmówców, że:

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. 

Ten mówi mi, że przed moim powyższym WriteArray (somearray) Mam ustawić zarówno na GeoTransform i SetProjection i SetUnitType na rzeczy do pracy (potencjalnie)

z samouczka GDAL API: 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 

moje przeprosiny za stworzenie zbyt długi post i spowiedź. Jeśli i kiedy dobrze to zrobię, dobrze byłoby mieć wszystkie notatki w jednym miejscu (długi post).

The Confession:

już wcześniej zrobione zdjęcie (JPEG), przekształcono go w GeoTIFF i importowane jako płytki do bazy danych PostGIS. Teraz próbuję utworzyć raster wysokościowy, na którym będzie zasłonięty obraz. To chyba wydaje się głupie, ale jest artysta, którego pragnę uczcić, podczas gdy na w tym samym czasie nie obrażam tych, którzy pracowali wytrwale nad stworzeniem i ulepszeniem tych wspaniałych narzędzi.

Artysta jest Belgiem, więc liczniki będą w porządku. Pracuje na dolnym Manhattanie w stanie Nowy Jork, UTM 18. Teraz rozsądna SetGeoTransform. Powyższe zdjęcie ma w = 3649/h = 2736, muszę się nad tym zastanowić.

Kolejna próba:

>>> 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 

Najwyraźniej coraz bliżej jednak jasne, czy SetUTM (18,1) została podchwycona. Czy to jest 4x4 w rzece Hudson lub Local_CS (układ współrzędnych) ? Co to jest cicha porażka?

więcej czytania

// 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. 

metrów 4x4 jest dość mała rozpiętość logiczne.

Być może jest to tak dobre, jak to tylko możliwe. SetGeoTransform pobiera jednostki w prawo, ustawia skalę i posiadasz przestrzeń Terragen World Space.

Ostatnia myśl: nie programuję, ale w pewnym zakresie mogę śledzić. Powiedział, że jeśli najpierw trochę zastanawiał, a potem co dane wyglądało w moim małym Terragen światowej przestrzeni (następujących dzięki http://www.gis.usu.edu/~chrisg/python/2009/ tydzień 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) 
>>> 

Więc to jest satysfakcjonujące. Wyobrażam sobie różnicę pomiędzy powyższym używanym numpy c , a ten wynik odnosi się do działań związanych z zastosowaniem Float16 w tym bardzo małym zakresie logicznym .

A na '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 
>>> 

Prawie całkowicie satysfakcjonujące, choć wygląda na to, że jestem w La Concordia Peru. Więc mam , aby dowiedzieć się, jak powiedzieć - jak bardziej na północ, jak New York North. Czy SetUTM przyjmuje trzeci element sugerujący "jak daleko" na północ lub na południe. Wygląda na to, że wczoraj natknąłem się na wykres UTM, który miał strefy etykiet listowych, coś jak C na równiku i może T lub S na obszarze Nowego Jorku.

Właściwie myślałem, że SetGeoTransform zasadniczo ustanawia lewy górny lewy i wschodni i to miało wpływ na to, jak daleko północ/południe "część SetUTM. Off to gdal-dev.

A potem jeszcze:

Miś Paddington udał się do Peru, bo miał bilet. Dotarłem tam, bo powiedziałem, że tam właśnie chciałem być. Terragen, działając w ten sposób, dał mi moją przestrzeń piksela. Kolejne wywołania srs działały na hf2 (SetUTM), ale wschód i północ zostały ustalone w Terragen, więc UTM 18 został ustawiony, ale w ramce ograniczającej na równiku. Wystarczająco dobry.

gdal_translate zabrał mnie do Nowego Jorku. Jestem w Windows, więc linia poleceń.i wynik;

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> 

Tak, wróciliśmy do NY. Prawdopodobnie istnieją lepsze sposoby podejścia do tego wszystkiego. I musiał mieć cel, przyjęty Tworzenie jak byłem postulując/improwizacji zestawu danych z numpy również. Muszę spojrzeć na inne formaty, które umożliwiają tworzenie. Elewacja w GeoTiff możliwość (myślę).

Moje podziękowania dla wszystkich uwag, sugestii i delikatnie trąca wobec odpowiedniego odczytu. Tworzenie map w Pythonie jest fajne!

Chris

Answer 1

nie są zbyt daleko. Twoim jedynym problemem są problemy z składnią w opcjach tworzenia GDAL. Wymienić:

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

z bez spacji przed lub po par klucz/wartość:

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

Sprawdź type(dst_ds) i powinno być <class 'osgeo.gdal.Dataset'> zamiast <type 'NoneType'> za cichy nie jak miało to miejsce w przypadku powyżej.

---

Aktualizacja Domyślnie warnings or exceptions are not shown. można je włączyć poprzez gdal.UseExceptions() zobaczyć co tyka pod spodem, np .:

>>> 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'>