创建海拔/高度字段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 import gdal 来自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文档的阅读器 作者：Daylon Graphics Ltd.的Ray Gardener * 该模块的部分由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 =                 （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教程中： 蟒蛇  导入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>
所以，我们回到了纽约。解决所有这些问题可能有更好的方法。一世 还必须有一个接受创建的目标，因为我也在从numpy推定/改进数据集。我需要查看允许创建的其他格式。我认为，GeoTiff中的高程是有可能的 感谢所有的评论，建议和对适当阅读的温柔推动。用python制作地图很有趣！ 克里斯     
2020-01-31 2 条评论
分享

没有找到相关结果

    已邀请:

    疼嘶桐

    您距离不太远。唯一的问题是GDAL创建选项的语法问题。更换： 
    [\'MINUSERPIXELVALUE = 1\',\'MAXUSERPIXELVALUE= 4\']
    键/值对之前或之后没有空格： 
    [\'MINUSERPIXELVALUE=1\', \'MAXUSERPIXELVALUE=4\']
    检查
    type(dst_ds)
    ，应该是
    <class \'osgeo.gdal.Dataset\'>
    ，而不是
    <type \'NoneType\'>
    。 更新默认情况下，不显示警告或异常。您可以通过ѭ14启用它们，以查看下面的滴答声，例如： 
    >>> 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\'>
        
    2020-01-31 0 0
    分享
    为什么被折叠? 0 个回复被折叠

    要回复问题请先登录或注册

    或代码 OrCode.com 备案号：粤ICP备15020848号-1
    Escape time: 0.075302839279175