创建高程/高度字段 gdal numpy python

Question

我想使用 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 驱动程序
  • 用途：Terragen TER 文档的阅读器
  • 作者：Ray Gardener，Daylon Graphics Ltd. *
  • 此模块的部分内容源自 GDAL 驱动程序
  • Frank Warmerdam，请参阅 http://www.gdal.org

我的歉意提前给雷·加德纳和弗兰克·沃默丹。

Terragen 格式注释：

用于写作： SCAL = 网格柱距离（米） hv_px = hv_m / SCAL span_px = span_m / SCAL 偏移量 = 参见 TerragenDataset::write_header() 比例 = 参见 TerragenDataset::write_header() 物理高压= (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)，将其转换为 geotiff，并将其作为图块导入到 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 世界空间。

最后的想法：我不编程，但在某种程度上我可以遵循。也就是说，我首先有点想知道我的小 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
>>> 

几乎完全令人满意，尽管看起来我在秘鲁的 La Concordia。所以我有 弄清楚怎么说——像更北的地方，比如纽约北。 SetUTM 是否采用第三个元素来表明向北或向南“多远”。昨天我好像看到了一张 UTM 图表，上面有字母标签区域，比如赤道上的 C，纽约地区可能是 T 或 S。

我实际上认为 SetGeoTransform 本质上是建立左上角的北向和东向，这影响了 SetUTM 的“北/南多远”部分。转到 gdal-dev。

再后来：

帕丁顿熊去了秘鲁，因为他有票。我到达那里是因为我说那就是我想去的地方。 Terragen 的工作方式为我提供了像素空间。随后对 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 中的提升是可能的（我认为）。

我感谢所有评论、建议以及对适当阅读的温和推动。用 python 制作地图很有趣！

克里斯

Answer 1

你离这个还不算太远。您唯一的问题是 GDAL 创建选项的语法问题。替换：

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

在键/值对之前或之后无空格：

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

检查type(dst_ds)，它应该是 而不是 来实现静默失败，如上述情况。

---

更新默认情况下，不显示警告或异常 。您可以通过 gdal.UseExceptions() 启用它们以查看下面的内容，例如：

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