Я хотел бы создать несколько растров высот / полей высоты, используя 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
Я полагаю, что упускаю что-то простое и жду ваших советов.
Спасибо
Chris
(продолжение позже)
- terragendataset.cpp, v 1.2
*
- Проект: 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
смещение = см. TerragenDataset :: write_header ()
масштаб = см. 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 для работы (потенциально)
Из учебника по GDAL API:
питон
импорт оср
импорт 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). Это 4x4 в
Река Гудзон или 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.
4x4 метра - довольно маленький логический промежуток.
Так что, возможно, это так же хорошо, как и получается. SetGeoTransform правильно определяет единицы измерения, устанавливает масштаб, и у вас есть Terragen World Space.
Заключительная мысль: я не программирую, но до такой степени, что могу следовать за ним. Тем не менее, сначала я подумал: а как эти данные выглядят в моем маленьком космическом мире Terragen?
(следующее спасибо 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
>>>
Почти полностью радует, хотя, похоже, я в Ла Конкордия Перу. Так что я
чтобы понять, как сказать - как север, как север Нью-Йорка. Берет ли SetUTM третий элемент, указывающий «как далеко» север или юг. Кажется, вчера я наткнулся на диаграмму UTM, в которой были зоны с буквенными метками, что-то вроде C на экваторе и, возможно, T или S для района Нью-Йорка.
Я на самом деле думал, что SetGeoTransform, по сути, устанавливал левый верхний север и восток, и это влияло на то, насколько далеко север / юг часть SetUTM. Выкл. В gdal-dev.
И еще позже:
Медведь Паддингтон отправился в Перу, потому что у него был билет.Я попал туда, потому что сказал, что хочу быть там.Терраген, работая как есть, дал мне место в пикселях.Последующие вызовы srs действовали на hf2 (SetUTM), но восток и север были установлены под Terragen, так что UTM 18 был установлен, но в ограничительной рамке на экваторе.Достаточно хорошо.
gdal_translate доставил меня в Нью-Йорк.Я в Windows, поэтому командная строка.и результат;
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 возможно (я думаю.)
Я благодарю все Комментарии, предложения и мягкие толчки к правильному чтению.Делать карты на питоне - это весело!
Крис