マップサーバーを作成するときに、高地の地形プロファイルを作成する必要がありました。 デジタル標高モデルのデータとして、SRTM（NASA Shuttle Radar Topography Mission）を使用することにしました。
代替のデータセットがありますが、SRTMは最も一般的であり、精度は非常に満足のいくものです。
データはいくつかのバージョンで提供されます-セルサイズが1秒角と3秒角のグリッド。 米国ではより正確な1秒データ（SRTM1）が利用できますが、地球の残りの部分（北緯60度から南緯54度の間）では3秒のデータ（SRTM3）しか利用できません。 ファイルは、1201x1201（または1秒バージョンの場合は3601x3601）値のマトリックスです。
公式ソースからデータを簡単にダウンロードするために、小さなスクリプト
をスケッチしました
srtm_download＃！/ bin / bash
使用法（）{
echo "使用法：" $ 0 "-region -lat_semisphere -lon_semisphere -min_lat -max_lat -min_lon -max_lon"
echo "地域の1つ：アフリカオーストラリアユーラシア諸島North_America South_America"
echo "例：srtm_download Eurasia NE 43 53 21 142"
1番出口
}
if（（ "$＃" <7））; それから
使い方
他に
地域= 1ドル
LATITUDE_SEMISPHERE = 2ドル
LONGITUDE_SEMISPHERE = 3ドル
MIN_LATITUDE = 4ドル
MAX_LATITUDE = 5ドル
MIN_LONGITUDE = 6ドル
MAX_LONGITUDE = 7ドル
SERVER_ADDRESS = "
dds.cr.usgs.gov/srtm/version2_1/SRTM3 $ {REGION}"
mkdir srtm
for i in $（seq $ MIN_LATITUDE $ MAX_LATITUDE）
する
$のj（seq $ MIN_LONGITUDE $ MAX_LONGITUDE）
する
wget -P srtm -nv "$ SERVER_ADDRESS / $ LATITUDE_SEMISPHERE $ i $ LONGITUDE_SEMISPHERE`printf％03d $ j`" "。hgt.zip"
やった
やった
fi
使用可能なすべてのDEMファイルは、srtmフォルダーに保存されます。
高地のプロファイルを作成するとき
、オランダのosmユーザー
Sjors Provoostの プロジェクトが基礎として
採用されました。
osm-serverタイルをインストールすることから始めました。 Ubuntu 12.04では、必要なすべてのコンポーネントを手動で組み立てなく
ても便利で迅速な
方法があります。 Debian向けのビルドでも、特別な問題は発生しませんでした。 これを行う方法は、
ここと
ここで見ることができ
ます 。
PostgreSQLは、これらのセットを保存および操作するために使用されます。 まず、PostgreSQLでsrtmというユーザーとデータベースを作成する必要があります
sudo -u postgres -icreateuser username # answer yes for superuser (although this isn't strictly necessary)createdb -E UTF8 -O username srtmexitpsycopg2（Python用のPostgreSQLアダプター）を使用してテーブルを作成してデータを入力し、GDALライブラリを使用してSRTMを読み取ると便利です。 以前に作成され、データで満たされたSRTMフォルダーで小さなPythonスクリプトを実行する必要があります。
srtm2pgsql＃！/ usr / bin / python
pgのインポート、psycopg2
osgeo import gdal、gdal_arrayから
輸入OS
zipファイルをインポート
再インポート
インポートシステム
数学インポートsqrtから
def getLatLonFromFileName（名前）：
＃latとlonに分割：
p = re.compile（ '[NSEW] \ d *'）
[lat_str、lon_str] = p.findall（名前）
＃北か南か？
lat_str [0] ==“ N”の場合：
lat = int（lat_str [1：]）
その他：
lat = -int（lat_str [1：]）
＃東か西か？
lon_str [0] ==“ E”の場合：
lon = int（lon_str [1：]）
その他：
lon = -int（lon_str [1：]）
return [lat、lon]
def loadTile（ファイル名）：
＃解凍する
zf = zipfile.ZipFile（ファイル名+ ".hgt.zip"）
zf.namelist（）の名前の場合：
outfile = open（名前、「wb」）
outfile.write（zf.read（名前））
outfile.flush（）
outfile.close（）
＃読む
srtm = gdal.Open（ファイル名+ '.hgt'）
＃クリーンアップ
os.remove（ファイル名+ '.hgt'）
return gdal_array.DatasetReadAsArray（srtm）
def posFromLatLon（lat、lon）：
return（lat * 360 + lon）* 1200 * 1200
クラスDatabasePg：
def __init __（self、db_name、db_user、db_pass）：
self.db = pg.DB（dbname = db_name、host = 'localhost'、user = db_user、passwd = db_pass）
def createTableAltitude（self）：
テーブル= self.db.get_tables（）
そうでない場合（表の「public.altitude」）：
self.db.query（ "\
CREATE TABLE高度（\
pos bigint NOT NULL、\
alt int NULL、\
主キー（pos）\
）; \
"）
真を返す
DatabasePsycopg2クラス：
def __init __（self、db_name、db_user、db_pass）：
self.db_name = db_name
conn = psycopg2.connect（ "dbname = '" + db_name + "' host = 'localhost' user = '" + db_user + "' password = '" + db_pass + "'"）
self.cursor = conn.cursor（）
def insertTile（self、tile、lat0、lon0）：
＃Psycopg2接続を使用し、そのcopy_toおよび
＃copy_fromコマンド。より効率的なCOPYコマンドを使用します。
＃この方法には一時ファイルが必要です。
＃開始位置を計算
begin = posFromLatLon（lat0、lon0）
＃最初に、データを一時ファイルに書き込みます。
f =オープン（ '/ tmp / tempcopy'、 'w'）
＃一番上の行と右の列をドロップします。
範囲内の行（1、len（タイル））：
範囲（0、len（タイル）-1）のcolの場合：
f.write（str（\
begin +（row-1）* 1200 + col \
）+ "\ t" + str（タイル[行] [列]）+ "\ n"）
f.close（）
＃次に、一時ファイルからデータを読み取り、
＃高度テーブル。
f =オープン（ '/ tmp / tempcopy'、 'r'）
サブプロセスのインポート
psqlcmd = "/ usr / bin / psql"
p = subprocess.Popen（ 'psql' + self.db_name + '-c“ STDINからの高度のコピー;”'、stdin = f、shell = True）;
p.wait（）
f.close
def main（）：
if len（sys.argv）！= 2：
「エラー！ \ n使用法：strm2pgsqlユーザー名»
sys.exit（）
その他：
db_pg = DatabasePg（ 'srtm'、sys.argv [1]、 ''）
db_psycopg2 = DatabasePsycopg2（ 'srtm'、sys.argv [1]、 ''）
db_pg.createTableAltitude（）
os.listdir内のファイル（ '。'）：
if（re.search（ 'hgt'、file））：
tilename = file.split（ '。'）[0]
「タイルの読み込み」+タイル名+「srtmデータベースへ」を印刷
[lat、lon] = getLatLonFromFileName（タイル名）
tile = loadTile（タイル名）
db_psycopg2.insertTile（タイル、緯度、経度）
__name__ == '__main__'の場合：
メイン（）
Qtを使用して直接高度プロファイルを作成します。 データベースとの接続を作成します（データベース
にアクセスするには
ドライバーが必要な場合があり
ます ）
db = QSqlDatabase::addDatabase("QPSQL"); db.setHostName("localhost"); db.setDatabaseName("srtm"); db.setUserName("user"); db.setPassword("pass"); bool ok = db.open();
標高プロファイルを作成するルートは、緯度、経度、高度を保存するGeoPointポイントで構成されます
geopoint.hクラスGeoPoint
{
公開：
GeoPoint（）;
GeoPoint（二重緯度、二重経度、二重高度）;
二重緯度（）;
二重経度（）;
二重高度（）;
void setLatitude（二重緯度）;
void setLongitude（二重経度）;
void setAltitude（二重高度）;
プライベート：
double m_latitude;
double m_longitude;
double m_altitude;
};
必要なポイント数を取得するには、ルート上の各ポイントの高さを補間して推定値を取得する必要があります。
interpolateRoute(route, numberRoutePoint); foreach (GeoPoint point, route) { point.setAltitude(altitude(db, point.latitude(), point.longitude())); }
補間のために、ルートの全長を見つけ（
Vincenti式を使用）、必要な数の新しいポイントをルートに追加します
補間 void interpolateRoute(QVector<GeoPoint> &route, int number) {
各ルートポイントの高さの推定値を取得するために、データベースで4つの最も近いポイントを見つけ、
双線形補間を実行します
高さ推定 double altitude(PostgreSqlDatabase db, double latitude, double longitude) { qulonglong tl, tr, bl, br; double a, b; posFromLatLon(latitude, longitude, tl, tr, bl, br, a, b); int atl = db.fetchAltitude(tl); int atr = db.fetchAltitude(tr); int abl = db.fetchAltitude(bl); int abr = db.fetchAltitude(br); return bilinearInterpolation(atl, atr, abl, abr, a, b); } int PostgreSqlDatabase::fetchAltitude(qulonglong pos) { QSqlQuery query; if (query.exec("SELECT * FROM altitude WHERE pos = " + QString::number(pos))) { while (query.next()) { return query.value(1).toInt(); } } return Geo::undefinedAltitude; } void posFromLatLon(const double &latitude, const double &longitude, qulonglong &tl, qulonglong &tr, qulonglong &bl, qulonglong &br, double &a, double &b) { double lat0 = floor(latitude); double lon0 = floor(longitude); qulonglong pos0 = (lat0*360 + lon0) * 1200*1200; double latPos = (1199./1200 - (latitude - lat0)) * 1200*1200; double lonPos = (longitude - lon0) * 1200; double latPosTop = floor(latPos / 1200) * 1200; double latPosBottom = ceil(latPos / 1200) * 1200; double lonPosLeft = floor(lonPos); double lonPosRight = ceil(lonPos); a = (latPos - latPosTop) / 1200; b = (lonPos - lonPosLeft); tl = qulonglong(pos0 + latPosTop + lonPosLeft); tr = qulonglong(pos0 + latPosTop + lonPosRight); bl = qulonglong(pos0 + latPosBottom + lonPosLeft); br = qulonglong(pos0 + latPosBottom + lonPosRight); } double bilinearInterpolation(const double &tl, const double &tr, const double &bl, const double &br, const double &a, const double &b) { double b1 = tl; double b2 = bl - tl; double b3 = tr - tl; double b4 = tl - bl - tr + br; return b1 + b2*a + b3*b + b4*a*b; }
したがって、出口には、ルートポイントのセットとその高さの推定値があります。
STRMデータの使用に関するトピックがコミュニティにとって興味深いものである場合、次の記事で、osmデータで使用するための色付きレリーフ基板の作成方法を共有できます。