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

从osgeo导入gdal
从gdalconst导入*

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

提前对Ray Gardener和Frank Warmerdam表示歉意。

Terragen格式说明：

写入时：
SCAL = 网格点距离（米）
hv_px = hv_m / SCAL
span_px = span_m / SCAL
offset = 见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
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 

抱歉发了这么长的帖子，还有一个小 confession。如果我能搞定这一切，把所有笔记放在一个地方（这个长帖子）会很不错。

小 confession：

我之前拍过一张照片（jpeg），把它转换成了geotiff，并作为切片导入到PostGIS数据库中。现在我想创建一个高程栅格图，以便在其上覆盖这张照片。这可能看起来有点傻，但我想向一位艺术家致敬，同时也不想冒犯那些努力创建和改进这些伟大工具的人。

这位艺术家是比利时人，所以用米作为单位是合适的。她在纽约下曼哈顿工作，所以使用UTM 18。现在需要设置一个合理的SetGeoTransform。上面提到的照片宽3649，高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（坐标系统）？什么是silent fail？

更多阅读

// 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/ 第四周）：

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

之后：

帕丁顿熊去秘鲁是因为他有票。我去那是因为我说我想去那里。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)，它应该显示为<class 'osgeo.gdal.Dataset'>，而不是<type 'NoneType'>，后者表示静默失败，就像上面那样。

---

更新 默认情况下，警告或异常不会显示。你可以通过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'>