標高/高さフィールドGDAL Numpy Pythonの作成
https://stackoverflow.com/questions/6319889
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
Russian質問
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からGDALCONSTインポートからインポートGDALから *
>>> 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
レイ・ガーデナーとフランク・ウォーマーダムに事前に謝罪します。
Terragen形式のメモ：
書き込みの場合：スカル=メーターのグリッドポスト距離hv_px = hv_m / scal span_px = span_m / scal offset = terragendataset :: write_header（）scale = terragendataset :: write_header（）物理hv =（hv_px -offset） * 6555556.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
過度に長い投稿と告白を作成したことをお詫びします。私がこれを正しくしたら、すべてのメモを1つの場所（長い投稿）に持っているといいでしょう。
告白：
私は以前に写真（JPEG）を撮り、それをジオティフに変換し、タイルとしてポストGISデータベースにインポートしました。私は今、写真をドレープするための標高ラスターを作成しようとしています。これはおそらくばかげているように思えますが、私が尊敬したいアーティストがいますが、同時にこれらの素晴らしいツールを作成し、改善するために熱心に働いた人々を怒らせません。
アーティストはベルギー人なので、メーターは整っています。彼女はニューヨーク州下部マンハッタンで働いています。上記の写真は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はユニットを正しく取得し、スケールを設定すると、テラゲンの世界空間があります。
最終的な考え：私はプログラムしていませんが、ある程度続くことができます。そうは言っても、私は最初に私の小さなテラゲンの世界空間でどのように見えるのか、そしてそれから私はちょっと疑問に思いました（次のことに感謝します 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は、北または南の「どこまで」を示唆する3番目の要素を取ります。昨日、EquatorのCのようなレターラベルゾーン、おそらくニューヨークエリアのTまたはSのようなものがあるUTMチャートに出くわしたようです。
私は実際に、SetGeoTransformが本質的に左上のnorthingを確立し、Eastingを確立していると思っていましたが、これはSetutmの「北/南/南」部分に影響を与えていました。 GDAL-DEVにオフ。
そして後で：
パディントンベアはチケットを持っていたのでペルーに行きました。私がそこに着いたのは、それが私がなりたかった場所だと言ったからです。テラゲンは、それがそうであるように働いて、私にピクセルのスペースを与えてくれました。その後のSRSへの呼び出しはHF2（setUtm）に作用しましたが、出掘とノーはテラゲンの下で確立されたため、UTM 18は設定されましたが、赤道の境界ボックスにありました。十分です。
GDAL_TRANSLATEは私をニューヨークに連れて行ってくれました。私は窓にいるのでコマンドラインです。そして結果;
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からデータセットを仮定/即興で見ているため、Createを受け入れたターゲットを持たなければなりませんでした。作成を許可する他の形式を見る必要があります。 Geotiffの標高は可能性があります（私は思う）
すべてのコメント、提案、適切な読書に対する穏やかな微調整に感謝します。 Pythonでマップを作るのは楽しいです！
クリス
解決
あなたはそれほど遠くない。唯一の問題は、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'>
