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
osgeoimportgdalからgdalconstimportから*
>>> 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ドライバー
- 目的：TerragenTERドキュメントのリーダー
- 著者：Ray Gardener、Daylon Graphics Ltd. *
- GDALドライバーから派生したこのモジュールの一部
- フランク・ウォーマーダム、http： //www.gdal.orgを参照
RayGardenerとFrankWarmerdamに事前にお詫び申し上げます。
Terragenフォーマットノート：
書き込みの場合：SCAL =グリッドポストの距離（メートル単位）hv_px = hv_m / SCAL span_px = span_m/SCALオフセット=TerragenDataset:: write_header（）を参照scale = 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
過度に長い投稿と告白を作成してしまったことをお詫びします。これが正しければ、すべてのメモを1か所（長い投稿）にまとめておくと便利です。
告白：
以前に写真（jpeg）を撮り、それをgeotiffに変換し、タイルとしてPostGISデータベースにインポートしました。現在、画像をドレープするための標高ラスターを作成しようとしています。これはおそらくばかげているように思えますが、私が尊敬したいアーティストがいますが、同時にこれらの優れたツールの作成と改善に熱心に取り組んできた人々を怒らせません。
アーティストはベルギー人なので、メーターが適切です。彼女はニューヨークのロウアーマンハッタン、UTM18で働いています。上記の写真は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 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)
>>>
ですから、これは満足のいくものです。上記で使用したnumpycとこの結果の違いは、この非常に小さな論理スパン全体に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番目の要素を取りますか。昨日、赤道でC、ニューヨークエリアでTまたはSのような文字ラベルゾーンがあるUTMチャートに出くわしたようです。
私は実際、SetGeoTransformが本質的に左上の北と東を確立していると思っていました。これは、SetUTMの「南北の距離」の部分に影響を与えていました。gdal-devに移動します。
そして後でまだ：
パディントンベアはチケットを持っていたのでペルーに行きました。そこに行きたいと言ったのでそこに着きました。Terragenは、そのように機能して、ピクセルスペースを与えてくれました。その後のsrsの呼び出しは、hf2（SetUTM）に作用しましたが、東と北はTerragenの下で確立されたため、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でマップを作成するのは楽しいです！
クリス