🤞🏼 👨🏽‍🏭 👩🏼 SRTMデータを使用したデジタル標高モデル 🚩 🛣️ 👨🏾‍🚒

マップサーバーを作成するときに、高地の地形プロファイルを作成する必要がありました。デジタル標高モデルのデータとして、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 -i

createuser username # answer yes for superuser (although this isn't strictly necessary)

createdb -E UTF8 -O username srtm

exit

psycopg2（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) { //     double length = routeLength(route); double minRes = length / number; //        : int i = 0; while (i < route.size() - 1) { QPair<GeoPoint, GeoPoint> pair = qMakePair(route[i], route[i+1]); int number = numberExtraPoints(pair, minRes); addPointsToRoute(route, i, number); i = i + number + 1; } } double routeLength(const QVector<GeoPoint> &route) { double length = 0; Geo geo; for (int i=0; i<route.size() - 1; i++) { GeoPoint origin = route[i]; GeoPoint destination = route[i+1]; length += geo.distance(origin, destination); } return length; } int numberExtraPoints(QPair<GeoPoint, GeoPoint> pair, double minRes) { Geo geo; double length = geo.distance(pair.first, pair.second); return qMax(int(floor(length / minRes) - 1), 0); } void addPointsToRoute(QVector<GeoPoint> &route, int index, int number) { double lat1 = route[index].latitude(); double lat2 = route[index+1].latitude(); double lon1 = route[index].longitude(); double lon2 = route[index+1].longitude(); double latRes = (lat2 - lat1) / number; double lonRes = (lon2 - lon1) / number; for (int i=0; i<number; i++) { lat1 += latRes; lon1 += lonRes; GeoPoint point(lat1, lon1, 0); route.insert(index + i, point); } } float WGS84MajorAxis = 6378.137; // km float WGS84MinorAxis = 6356.7523142; // km float WGS84Flattering = 1 / 298.257223563; double distance(GeoPoint origin, GeoPoint destination) { double lat1 = origin.latitude(); double lat2 = destination.latitude(); double lng1 = origin.longitude(); double lng2 = destination.longitude(); double major = WGS84MajorAxis; double minor = WGS84MinorAxis; double f = WGS84Flattering; double deltaLng = lng2 - lng1; double reducedLat1 = atan((1 - f) * tan(lat1)); double reducedLat2 = atan((1 - f) * tan(lat2)); double sinReduced1 = sin(reducedLat1); double sinReduced2 = sin(reducedLat2); double cosReduced1 = cos(reducedLat1); double cosReduced2 = cos(reducedLat2); double lambdaLng = deltaLng; double lambdaPrime = 2 * M_PI; double iterLimit = 20; double sinSigma = 0; double cosSigma = 0; double cos2SigmaM = 0; double sigma = 0; double cosSqAlpha = 0; while ((abs(lambdaLng - lambdaPrime) > 10e-12) && (iterLimit > 0)) { double sinLambdaLng = sin(lambdaLng); double cosLambdaLng = cos(lambdaLng); sinSigma = sqrt((cosReduced2*sinLambdaLng) * (cosReduced2*sinLambdaLng) + (cosReduced1*sinReduced2 - sinReduced1*cosReduced2*cosLambdaLng) * (cosReduced1*sinReduced2 - sinReduced1*cosReduced2*cosLambdaLng)); if (sinSigma == 0) return 0; //   cosSigma = sinReduced1*sinReduced2 + cosReduced1*cosReduced2*cosLambdaLng; sigma = atan2(sinSigma, cosSigma); double sinAlpha = cosReduced1*cosReduced2*sinLambdaLng/sinSigma; cosSqAlpha = 1 - sinAlpha*sinAlpha; cos2SigmaM = 0.; //   if (cosSqAlpha != 0) { cos2SigmaM = cosSigma - 2*(sinReduced1*sinReduced2/cosSqAlpha); } double C = f / 16. * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha)); lambdaPrime = lambdaLng; lambdaLng = (deltaLng + (1 - C)*f*sinAlpha*(sigma + C * sinSigma * (cos2SigmaM + C * cosSigma *(-1 + 2*cos2SigmaM*cos2SigmaM)))); --iterLimit; } if (iterLimit == 0) { return 0; log->info() << "Geo, calculate distance: Vincenty formula failed to converge!"; } double uSq = cosSqAlpha * (major*major - minor*minor) / (minor*minor); double A = 1 + uSq / 16384. * (4096 + uSq * (-768 + uSq *(320 - 175 * uSq))); double B = uSq / 1024. * (256 + uSq * (-128 + uSq * (74 - 47 * uSq))); double deltaSigma = (B * sinSigma *(cos2SigmaM + B / 4. *(cosSigma * (-1 + 2 * cos2SigmaM*cos2SigmaM) - B / 6. * cos2SigmaM * (-3 + 4 * sinSigma*sinSigma)*(-3 + 4 * cos2SigmaM*cos2SigmaM)))); return minor * A * (sigma - deltaSigma); }

各ルートポイントの高さの推定値を取得するために、データベースで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データで使用するための色付きレリーフ基板の作成方法を共有できます。

SRTMデータを使用したデジタル標高モデル

More articles: