SRTM डेटा का उपयोग करके डिजिटल उन्नयन मॉडल

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