создание поля высоты / высоты gdal numpy python

Я хотел бы создать несколько растров высот / полей высот, используя python, gdal и numpy. Я застрял на numpy (и, вероятно, на python и gdal.)

В numpy я пытался сделать следующее:

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

из osgeo import gdal из 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

Я полагаю, что мне не хватает чего-то простого, и с нетерпением жду вашего совета.

Спасибо,

Крис

(продолжение позже)

• terragendataset.cpp,v 1.2 *
  • Project: Terragen(tm) TER Driver
  • Назначение: Читатель документов Terragen TER
  • Автор: Рэй Гарденер, Daylon Graphics Ltd. *
  • Части этого модуля, полученные из драйверов GDAL,
  • Фрэнк Вармердам, см. http://www.gdal.org

Заранее приношу свои извинения Рэю Гарденеру и Фрэнку Вармердаму.

Примечания к формату Terragen:

Для записи: SCAL = расстояние по сетке в метрах hv_px = hv_m / SCAL span_px = span_m / SCAL offset = см. TerragenDataset :: write_header () scale = см. TerragenDataset :: write_header () физический hv = (hv_px - смещение) * 65536.0 / масштаб

Мы сообщаем вызывающим абонентам, что:

    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.

Это говорит мне, что до моего вышеописанного WriteArray (somearray) я должен установить как GeoTransform, так и SetProjection и SetUnitType, чтобы все работало (потенциально)

Из учебного пособия по API GDAL: 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 

Приношу свои извинения за слишком длинный пост и признание. Если и когда я сделаю это правильно, было бы неплохо собрать все заметки в одном месте (длинный пост).

Признание:

Я ранее сделал снимок (jpeg), преобразовал его в геотиф и импортировал его как плитки в базу данных PostGIS. Сейчас я пытаюсь создать растр высот, чтобы наложить на него изображение. Возможно, это кажется глупым, но есть художник, которого я хочу воздать должное, но в то же время не оскорбить тех, кто усердно трудился над созданием и улучшением этих великих инструментов.

Художник бельгиец, поэтому метры в порядке. Она работает в нижнем Манхэттене, штат Нью-Йорк, так, UTM 18. Теперь какой-нибудь разумный SetGeoTransform. Вышеупомянутое изображение имеет вид w = 3649 / h = 2736, я должен над этим подумать.

Еще одна попытка:

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

Явно приближаюсь, но неясно, был ли поднят SetUTM (18,1). Это внедорожник на реке Гудзон или Local_CS (система координат)? Что за тихий провал?

Больше Чтения

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

4х4 метра - довольно маленький логический промежуток.

Так что, возможно, это так хорошо, как есть. SetGeoTransform правильно определяет единицы измерения, устанавливает масштаб, и у вас есть ваше мировое пространство Terragen.

Заключительная мысль: я не программирую, но до некоторой степени могу следовать за мной. Тем не менее, я сначала подумал, а потом, как данные выглядели в моем маленьком Terragen World Space (следующее спасибо http://www.gis.usu.edu/~chrisg/python/2009/ неделя 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)
>>>

Так что это отрадно. Я предполагаю, что разница между использованным выше numpy c и этим результатом относится к действиям по применению Float16 в этом очень маленьком логическом интервале.

А теперь о 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
>>> 

Почти полностью отрадно, хотя похоже, что я нахожусь в La Concordia Peru. Так что я должен придумать, как сказать - вроде севернее, как север Нью-Йорка. Принимает ли SetUTM третий элемент, предлагающий «как далеко» к северу или югу. Кажется, вчера я наткнулся на диаграмму UTM, на которой были обозначены зоны букв, что-то вроде C на экваторе и, возможно, T или S для района Нью-Йорка.

На самом деле я думал, что SetGeoTransform по существу устанавливает верхний левый север и восток, и это влияло на ту часть SetUTM, «как далеко север / юг». Переходим на gdal-dev.

И еще позже:

Паддингтон Медведь отправился в Перу, потому что у него был билет. Я попал туда, потому что сказал, что хочу быть именно там. Terragen, как он работает, предоставил мне мое пиксельное пространство. Последующие вызовы srs действовали на hf2 (SetUTM), но восток и север были установлены под Terragen, поэтому UTM 18 был установлен, но в ограничивающей рамке на экваторе. Достаточно хорошо.

gdal_translate доставил меня в Нью-Йорк. Я в окнах так командная строка. и результат;

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>

Итак, мы вернулись в Нью-Йорк. Наверное, есть лучшие способы подойти ко всему этому. У меня должна была быть цель, которая принимала Create, поскольку я также постулировал / импровизировал набор данных из numpy. Мне нужно посмотреть другие форматы, позволяющие создавать. Подъем в GeoTiff возможен (я думаю).

Моя благодарность за все комментарии, предложения и легкие подталкивания к правильному чтению. Создавать карты на Python - это весело!

Крис