создание поля высоты / высоты 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

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

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

Спасибо,

Крис

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

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

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

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

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

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

    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, чтобы все работало (потенциально)

Из учебного руководства GDAL API: Python импорт 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)

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

>>> 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). Это 4х4 в Река Гудзон или 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 «как далеко север / юг». Поехали в гдал-дев.

И еще позже:

Медведь Паддингтон отправился в Перу, потому что у него был билет. Я попал туда, потому что сказал, что хочу быть именно там. 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 - это весело!

Крис