python - 创建高程/高度场 gdal numpy python

我想使用 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 - 偏移) * 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 教程:
Python
导入操作系统
导入 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 在
hudson 河还是 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 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。

后来还是:

帕丁顿熊去秘鲁是因为他有一张票。我到那里是因为我说过那是我想去的地方。 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制作 map 很有趣!

克里斯

最佳答案

你离得不远了。您唯一的问题是 GDAL 创建选项的语法问题。代替:

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

在键/值对之前或之后没有空格:

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

查询 type(dst_ds)应该是 <class 'osgeo.gdal.Dataset'>而不是 <type 'NoneType'>对于上述情况的无声失败。

更新默认情况下，warnings or exceptions are not shown .您可以通过 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'>

关于python - 创建高程/高度场 gdal numpy python，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/6319889/

python - 创建高程/高度场 gdal numpy python

上一篇：python - 如何在 wxpython 中调整我的工具提示？

下一篇：使用 pyPDF 删除空白页的 Python 脚本