Digital elevation model using SRTM data
When creating a map server, I needed to create high-altitude terrain profiles. I decided to use SRTM (NASA Shuttle Radar Topography Mission) as a data for a digital elevation model. Although there are alternative data sets , SRTM is the most common and quite satisfied me in accuracy. The data is supplied in several versions - a grid with a cell size of 1 arc-second and 3 arc-seconds. More accurate one-second data (SRTM1) is available on the territory of the United States, only three-second data (SRTM3) is available on the rest of the earth’s surface (between 60 degrees north latitude and 54 degrees south latitude). The files are a matrix of 1201x1201 (or 3601x3601 for the one-second version) values.
For convenient downloading of data from an official source, a small script is sketched.
All available files with DEM will be saved in the srtm folder.
The construction of high-altitude profiles was based on the project of the Dutch osm-user Sjors Provoost .
I started with the installation of an osm server tiles. Under Ubuntu 12.04 there is a convenient and fast way to do without manually assembling all the necessary components. Build under Debian, too, did not cause any problems. How to do this can be found here and here .
PostgreSQL is used to store and work with these sets. First we need to create a user in PostgreSQL and a database named srtm
')
It is convenient to use psycopg2 (PostgreSQL adapter for Python) to create and populate the table, and GDAL library for reading SRTM. A small python script that needs to be run in the SRTM folder that was previously created and filled with data.
I used Qt to create a high profile directly. Create a connection with the database (you may need a driver to access the database )
The route whose high-altitude profile we want to create consists of GeoPoint points that store latitude, longitude and height
To obtain the required number of points, you need to interpolate and get a height estimate for each point of the route.
For interpolation I find the total length of the route (using the Vincenti's formula) and add the required number of new points to the route.
To obtain an estimate of the height of each point of the route, we find the four nearest points in the database and perform bilinear interpolation
Thus at the exit we will have a set of route points with an estimate of their height.
If the topic of using STRM data is interesting to the community, in the next article I can share with you how I made a color relief substrate for use with osm data.
For convenient downloading of data from an official source, a small script is sketched.
srtm_download
#! / bin / bash
usage () {
echo "Usage:" $ 0 "-region -lat_semisphere -lon_semisphere -min_lat -max_lat -min_lon -max_lon"
echo "Region one of: Africa Australia Eurasia Islands North_America South_America"
echo "Example: srtm_download Eurasia NE 43 53 21 142"
exit 1
}
if (("$ #" <7)); then
usage
else
REGION = $ 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)
do
for j in $ (seq $ MIN_LONGITUDE $ MAX_LONGITUDE)
do
wget -P srtm -nv "$ SERVER_ADDRESS / $ LATITUDE_SEMISPHERE $ i $ LONGITUDE_SEMISPHERE`printf% 03d $ j`" ". hgt.zip"
done
done
fi
usage () {
echo "Usage:" $ 0 "-region -lat_semisphere -lon_semisphere -min_lat -max_lat -min_lon -max_lon"
echo "Region one of: Africa Australia Eurasia Islands North_America South_America"
echo "Example: srtm_download Eurasia NE 43 53 21 142"
exit 1
}
if (("$ #" <7)); then
usage
else
REGION = $ 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)
do
for j in $ (seq $ MIN_LONGITUDE $ MAX_LONGITUDE)
do
wget -P srtm -nv "$ SERVER_ADDRESS / $ LATITUDE_SEMISPHERE $ i $ LONGITUDE_SEMISPHERE`printf% 03d $ j`" ". hgt.zip"
done
done
fi
All available files with DEM will be saved in the srtm folder.
The construction of high-altitude profiles was based on the project of the Dutch osm-user Sjors Provoost .
I started with the installation of an osm server tiles. Under Ubuntu 12.04 there is a convenient and fast way to do without manually assembling all the necessary components. Build under Debian, too, did not cause any problems. How to do this can be found here and here .
PostgreSQL is used to store and work with these sets. First we need to create a user in PostgreSQL and a database named srtm
')
sudo -u postgres -icreateuser username # answer yes for superuser (although this isn't strictly necessary)createdb -E UTF8 -O username srtmexitIt is convenient to use psycopg2 (PostgreSQL adapter for Python) to create and populate the table, and GDAL library for reading SRTM. A small python script that needs to be run in the SRTM folder that was previously created and filled with data.
srtm2pgsql
#! / usr / bin / python
import pg, psycopg2
from osgeo import gdal, gdal_array
import os
import zipfile
import re
import sys
from math import sqrt
def getLatLonFromFileName (name):
# Split up in lat and lon:
p = re.compile ('[NSEW] \ d *')
[lat_str, lon_str] = p.findall (name)
# North or south?
if lat_str [0] == "N":
lat = int (lat_str [1:])
else:
lat = -int (lat_str [1:])
# East or west?
if lon_str [0] == "E":
lon = int (lon_str [1:])
else:
lon = -int (lon_str [1:])
return [lat, lon]
def loadTile (filename):
# Unzip it
zf = zipfile.ZipFile (filename + ".hgt.zip")
for name in zf.namelist ():
outfile = open (name, 'wb')
outfile.write (zf.read (name))
outfile.flush ()
outfile.close ()
# Read it
srtm = gdal.Open (filename + '.hgt')
# Clean up
os.remove (filename + '.hgt')
return gdal_array.DatasetReadAsArray (srtm)
def posFromLatLon (lat, lon):
return (lat * 360 + lon) * 1200 * 1200
class 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):
tables = self.db.get_tables ()
if not ('public.altitude' in tables):
self.db.query ("\
CREATE TABLE altitude (\
pos bigint NOT NULL, \
alt int NULL, \
PRIMARY KEY (pos) \
); \
")
return true
class 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):
# I use the Psycopg2 connection, with its copy_to and
# copy_from commands, which use the more efficient COPY command.
# This method requires a temporary file.
# Calculate begin position
begin = posFromLatLon (lat0, lon0)
# First we write the data into a temporary file.
f = open ('/ tmp / tempcopy', 'w')
# We drop the top row and right column.
for row in range (1, len (tile)):
for col in range (0, len (tile) - 1):
f.write (str (\
begin + (row-1) * 1200 + col \
) + "\ t" + str (tile [row] [col]) + "\ n")
f.close ()
# Now we read the data file.
# altitude table.
f = open ('/ tmp / tempcopy', 'r')
import subprocess
psqlcmd = "/ usr / bin / psql"
p = subprocess.Popen ('psql' + self.db_name + '-c "COPY altitude FROM STDIN;"', stdin = f, shell = True);
p.wait ()
f.close
def main ():
if len (sys.argv)! = 2:
print “Error! \ nUsage: strm2pgsql username "
sys.exit ()
else:
db_pg = DatabasePg ('srtm', sys.argv [1], '')
db_psycopg2 = DatabasePsycopg2 ('srtm', sys.argv [1], '')
db_pg.createTableAltitude ()
for file in os.listdir ('.'):
if (re.search ('hgt', file)):
tilename = file.split ('.') [0]
print "Load tile" + tilename + "to srtm database"
[lat, lon] = getLatLonFromFileName (tilename)
tile = loadTile (tilename)
db_psycopg2.insertTile (tile, lat, lon)
if __name__ == '__main__':
main ()
import pg, psycopg2
from osgeo import gdal, gdal_array
import os
import zipfile
import re
import sys
from math import sqrt
def getLatLonFromFileName (name):
# Split up in lat and lon:
p = re.compile ('[NSEW] \ d *')
[lat_str, lon_str] = p.findall (name)
# North or south?
if lat_str [0] == "N":
lat = int (lat_str [1:])
else:
lat = -int (lat_str [1:])
# East or west?
if lon_str [0] == "E":
lon = int (lon_str [1:])
else:
lon = -int (lon_str [1:])
return [lat, lon]
def loadTile (filename):
# Unzip it
zf = zipfile.ZipFile (filename + ".hgt.zip")
for name in zf.namelist ():
outfile = open (name, 'wb')
outfile.write (zf.read (name))
outfile.flush ()
outfile.close ()
# Read it
srtm = gdal.Open (filename + '.hgt')
# Clean up
os.remove (filename + '.hgt')
return gdal_array.DatasetReadAsArray (srtm)
def posFromLatLon (lat, lon):
return (lat * 360 + lon) * 1200 * 1200
class 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):
tables = self.db.get_tables ()
if not ('public.altitude' in tables):
self.db.query ("\
CREATE TABLE altitude (\
pos bigint NOT NULL, \
alt int NULL, \
PRIMARY KEY (pos) \
); \
")
return true
class 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):
# I use the Psycopg2 connection, with its copy_to and
# copy_from commands, which use the more efficient COPY command.
# This method requires a temporary file.
# Calculate begin position
begin = posFromLatLon (lat0, lon0)
# First we write the data into a temporary file.
f = open ('/ tmp / tempcopy', 'w')
# We drop the top row and right column.
for row in range (1, len (tile)):
for col in range (0, len (tile) - 1):
f.write (str (\
begin + (row-1) * 1200 + col \
) + "\ t" + str (tile [row] [col]) + "\ n")
f.close ()
# Now we read the data file.
# altitude table.
f = open ('/ tmp / tempcopy', 'r')
import subprocess
psqlcmd = "/ usr / bin / psql"
p = subprocess.Popen ('psql' + self.db_name + '-c "COPY altitude FROM STDIN;"', stdin = f, shell = True);
p.wait ()
f.close
def main ():
if len (sys.argv)! = 2:
print “Error! \ nUsage: strm2pgsql username "
sys.exit ()
else:
db_pg = DatabasePg ('srtm', sys.argv [1], '')
db_psycopg2 = DatabasePsycopg2 ('srtm', sys.argv [1], '')
db_pg.createTableAltitude ()
for file in os.listdir ('.'):
if (re.search ('hgt', file)):
tilename = file.split ('.') [0]
print "Load tile" + tilename + "to srtm database"
[lat, lon] = getLatLonFromFileName (tilename)
tile = loadTile (tilename)
db_psycopg2.insertTile (tile, lat, lon)
if __name__ == '__main__':
main ()
I used Qt to create a high profile directly. Create a connection with the database (you may need a driver to access the database )
db = QSqlDatabase::addDatabase("QPSQL"); db.setHostName("localhost"); db.setDatabaseName("srtm"); db.setUserName("user"); db.setPassword("pass"); bool ok = db.open(); The route whose high-altitude profile we want to create consists of GeoPoint points that store latitude, longitude and height
geopoint.h
class GeoPoint
{
public:
GeoPoint ();
GeoPoint (double latitude, double longitude, double altitude);
double latitude ();
double longitude ();
double altitude ();
void setLatitude (double latitude);
void setLongitude (double longitude);
void setAltitude (double altitude);
private:
double m_latitude;
double m_longitude;
double m_altitude;
};
{
public:
GeoPoint ();
GeoPoint (double latitude, double longitude, double altitude);
double latitude ();
double longitude ();
double altitude ();
void setLatitude (double latitude);
void setLongitude (double longitude);
void setAltitude (double altitude);
private:
double m_latitude;
double m_longitude;
double m_altitude;
};
To obtain the required number of points, you need to interpolate and get a height estimate for each point of the route.
interpolateRoute(route, numberRoutePoint); foreach (GeoPoint point, route) { point.setAltitude(altitude(db, point.latitude(), point.longitude())); } For interpolation I find the total length of the route (using the Vincenti's formula) and add the required number of new points to the route.
Interpolation
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); } To obtain an estimate of the height of each point of the route, we find the four nearest points in the database and perform bilinear interpolation
Height estimate
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; } Thus at the exit we will have a set of route points with an estimate of their height.
If the topic of using STRM data is interesting to the community, in the next article I can share with you how I made a color relief substrate for use with osm data.
Source: https://habr.com/ru/post/161201/
All Articles