【发布时间】：2011-09-13 06:37:22
【问题描述】：
我想使用 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 驱动程序，由
- 弗兰克·沃默丹，见http://www.gdal.org
我提前向 Ray Gardener 和 Frank Warmerdam 道歉。
Terragen 格式说明：
写作： SCAL = 以米为单位的网格柱距离 hv_px = hv_m / SCAL span_px = span_m / SCAL 偏移量 = 见 TerragenDataset::write_header() scale = 参见 TerragenDataset::write_header() 物理hv = (hv_px - 偏移量) * 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 导入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制作地图很有趣！
克里斯
【问题讨论】：
-
问题是什么？您正在将零写入文件...您想做什么？写零不起作用吗？如果要将 numpy 数组写入文件，请将其传入，而不是您正在创建的 zeros 数组。 （您可能需要传入
c.astype(numpy.float32)，如果您正在创建一个 32 位波段并且c是一个 64 位数组（这取决于您使用的平台）。 -
The driver 只支持
gdal.GDT_Int16——不支持float32 -
@JoeKingston - 零值有效，但我确实想将 c 作为浮点数 32 传递，因为我正在创建自己的高度数据。