Цифровая модель рельефа с использованием SRTM данных

#!/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 N E 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

#!/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 from the temporary file and put it in the
# 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()

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;
};

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);
}

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;
}