我希望使用Python、GDAL和NumPy创建一些高程/高度栅格图。我在NumPy上遇到了困难（也可能是Python和GDAL）。
在NumPy中，我一直在尝试以下内容：
我觉得我错过了一些简单的东西，期待你的建议。
谢谢，
克里斯
（稍后继续）
terragendataset.cpp，v1.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 偏移量=参见TerragenDataset :: write_header() 比例尺=参见TerragenDataset :: write_header() 物理hv = （hv_px - offset）* 65536.0 / scale 我们告诉调用者：
抱歉我发了一篇过长的帖子和一个自白。如果我做对了，把所有笔记放在一个地方（长贴）会很好。
自白：
我之前拍了一张照片（jpeg），将其转换为geotiff，并将其作为瓦片导入到PostGIS数据库中。现在我正在尝试创建高程栅格以覆盖这张图片。这可能看起来很傻，但有位艺术家我想要表彰，同时又不冒犯那些努力创造和改进这些伟大工具的人们。
这位艺术家是比利时人，所以应该使用米为单位。她在纽约下曼哈顿，UTM 18。现在需要合理的SetGeoTransform。上述提到的图片是w=3649/h=2736，我需要考虑一下。
几乎完全令人满意，尽管看起来像我在秘鲁拉孔科迪亚。所以我必须想办法说-像更北，像纽约北部。SetUTM是否需要第三个元素来表示向北或向南的距离。昨天似乎遇到了一个包含字母标签区域的UTM图表，赤道上可能有C，纽约地区可能有T或S之类的东西。
实际上，我认为SetGeoTransform基本上是建立左上角的北向和东向，这影响了SetUTM中的“向北/向南多远”的部分。前往gdal-dev。
后来：
Paddington熊去了秘鲁，因为他有一张票。我到那里是因为我说那是我想去的地方。Terragen，按照它的工作方式，给了我我的像素空间。对srs的后续调用作用于hf2（SetUTM），但是东向和北向是在Terragen下建立的，因此UTM 18被设置在赤道上的边界框中。够好。
gdal_translate带我去了纽约。我在Windows上使用命令行。结果如下：
在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，v1.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 偏移量=参见TerragenDataset :: write_header() 比例尺=参见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
抱歉我发了一篇过长的帖子和一个自白。如果我做对了，把所有笔记放在一个地方（长贴）会很好。
自白：
我之前拍了一张照片（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还是本地坐标系？什么是无声失败？
更多阅读
// 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 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
>>>
几乎完全令人满意，尽管看起来像我在秘鲁拉孔科迪亚。所以我必须想办法说-像更北，像纽约北部。SetUTM是否需要第三个元素来表示向北或向南的距离。昨天似乎遇到了一个包含字母标签区域的UTM图表，赤道上可能有C，纽约地区可能有T或S之类的东西。
实际上，我认为SetGeoTransform基本上是建立左上角的北向和东向，这影响了SetUTM中的“向北/向南多远”的部分。前往gdal-dev。
后来：
Paddington熊去了秘鲁，因为他有一张票。我到那里是因为我说那是我想去的地方。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>
所以，我们回到了纽约。有可能有更好的方法来处理这一切。我必须要有一个目标，接受从numpy中推导/即兴创作的数据集，因此需要查看其他允许创建的格式。GeoTiff中的高程是一个可能性（我认为）。感谢所有的评论、建议和对适当阅读的温柔推动。用Python制作地图很有趣！
Chris
c.astype(numpy.float32))。 - Joe Kingtongdal.GDT_Int16，不支持 float32。 - Mike T