384 lines
11 KiB
C++
384 lines
11 KiB
C++
/*
|
|
* (c) 2002-2016 Hannes Krueger
|
|
* This file is part of the GPLIGC/ogie package
|
|
*
|
|
* This program is free software: you can redistribute it and/or modify
|
|
* it under the terms of the GNU General Public License as published by
|
|
* the Free Software Foundation, either version 3 of the License, or
|
|
* (at your option) any later version.
|
|
*
|
|
* This program is distributed in the hope that it will be useful,
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
* GNU General Public License for more details.
|
|
*
|
|
* You should have received a copy of the GNU General Public License
|
|
* along with this program. If not, see <http://www.gnu.org/licenses/>
|
|
*
|
|
*/
|
|
|
|
/*
|
|
This program merges a WORLD.DEM (created with createworld from the GTOPO30-tiles)
|
|
and an etopo2.i2 into a WORLD3.DEM. (world3 is the new version with water-flag!)
|
|
(WORLD2 was without waterflag)
|
|
All files have 2 byte integer-values in big-endian order.
|
|
|
|
Another option is to prepare subsets of GTOPO30 WORLDx.DEM file
|
|
--latmax --latmin --lonmax --lonmin --out file
|
|
|
|
*/
|
|
|
|
#include <fstream>
|
|
#include <iostream>
|
|
#include <cstring>
|
|
#include <cmath>
|
|
#include <cstdio>
|
|
#include <cstdlib>
|
|
|
|
//const double GTOPO30GRID = 0.00833333333333333;
|
|
//const double ETOPO2GRID = 0.03333333333333333;
|
|
const double GTOPO30GRID = 1.0/120.0;
|
|
const double ETOPO2GRID = 1.0/30.0;
|
|
|
|
using namespace std;
|
|
//const string rcsid_etopo2merger_cpp =
|
|
// "$Id: etopo2merger.cpp 3 2014-07-31 09:59:20Z kruegerh $";
|
|
|
|
short int etopo_array[5401][10801];
|
|
|
|
// This function uses bilinear interpolation to map values
|
|
// from etopo2 to gtopo30
|
|
short int getalt(double lat, double lon, int debug)
|
|
{
|
|
int steps_lat, steps_lon;
|
|
short int alt, h1, h2, h3, h4, max, min;
|
|
double lat_frac, lon_frac, tmp;
|
|
|
|
double dem_lon_min = -180.0;
|
|
double dem_lat_max = 90.0;
|
|
|
|
// Layout:
|
|
// h3 h4
|
|
// h1 h2
|
|
// lat_frac is the remaining fraction on the side h1 towards h3
|
|
// lon_frac is the remaining fraction on the side h1 towards h2
|
|
|
|
steps_lat = (int) (1.0 + (dem_lat_max - lat) / ETOPO2GRID);
|
|
steps_lon = (int) ((lon - dem_lon_min) / ETOPO2GRID);
|
|
|
|
lat_frac = -1.0 * (-steps_lat * ETOPO2GRID + dem_lat_max - lat) / ETOPO2GRID;
|
|
lon_frac = (-steps_lon * ETOPO2GRID - dem_lon_min + lon) / ETOPO2GRID;
|
|
h1 = etopo_array[steps_lat][steps_lon];
|
|
|
|
if(lat_frac < 0 || lon_frac < 0 || debug) printf("latfrac = %f, lonfrac = %f\n", lat_frac, lon_frac);
|
|
|
|
steps_lat = (int) (1.0 + (dem_lat_max - lat) / ETOPO2GRID);
|
|
steps_lon = (int) (1.0 + (lon - dem_lon_min) / ETOPO2GRID);
|
|
h2 = etopo_array[steps_lat][steps_lon];
|
|
|
|
steps_lat = (int) ((dem_lat_max - lat) / ETOPO2GRID);
|
|
steps_lon = (int) ((lon - dem_lon_min) / ETOPO2GRID);
|
|
h3 = etopo_array[steps_lat][steps_lon];
|
|
|
|
steps_lat = (int) ((dem_lat_max - lat) / ETOPO2GRID);
|
|
steps_lon = (int) (1.0 + (lon - dem_lon_min) / ETOPO2GRID);
|
|
h4 = etopo_array[steps_lat][steps_lon];
|
|
|
|
tmp = ((1.0 - lon_frac) * (1.0 - lat_frac) * (double) h1
|
|
+ lon_frac * (1.0 - lat_frac) * (double) h2
|
|
+ (1.0 - lon_frac) * lat_frac * (double) h3
|
|
+ lat_frac * lon_frac * (double) h4);
|
|
|
|
tmp = (double)((float) tmp); //This is in order to achieve identical results on x86_64 and i386
|
|
// Looks like Voodoo, but it can be understood: We are always sitting right on the fence
|
|
// rounding this way or that. If x86_64 works with different internal precision, then the
|
|
// rounding outcome in the cast operation (below) can go different ways...
|
|
|
|
if(tmp > 0.0){
|
|
alt = (short int) (0.5 + tmp);
|
|
} else {
|
|
alt = (short int) (-0.5 + tmp);
|
|
}
|
|
|
|
// Sanity checks
|
|
max = -10000;
|
|
min = 16000;
|
|
if(h1 > max) max = h1;
|
|
if(h2 > max) max = h2;
|
|
if(h3 > max) max = h3;
|
|
if(h4 > max) max = h4;
|
|
if(h1 < min) min = h1;
|
|
if(h2 < min) min = h2;
|
|
if(h3 < min) min = h3;
|
|
if(h4 < min) min = h4;
|
|
if(alt > max || alt < min || debug){
|
|
printf("Interpolation: h1 = %d, h2 = %d, h3 = %d, h4 = %d, alt = %d (%f), lat_frac = %f, lon_frac = %f\n",
|
|
h1, h2, h3, h4, alt, tmp, lat_frac, lon_frac);
|
|
}
|
|
// end sanity checks
|
|
|
|
return alt;
|
|
}
|
|
|
|
|
|
void failed(void)
|
|
{
|
|
cout << "Some of the input files are missing, "
|
|
<< "or we are unable to write to the output file."
|
|
<< endl
|
|
<<
|
|
"The files WORLD.DEM should be in this directory. And for etopo2merging etopo2.i2!"
|
|
<< endl;
|
|
exit(1);
|
|
}
|
|
|
|
//---------------------------------------------
|
|
|
|
int main(int argc, char *argv[])
|
|
{
|
|
bool ETOPOMERGING = true;
|
|
double minlat, maxlat, minlon, maxlon;
|
|
string filename;
|
|
string source = "WORLD.DEM";
|
|
|
|
// parsing of options...
|
|
for (int i = 0; i < argc; i++) {
|
|
|
|
if (strcmp(argv[i], "--latmin") == 0)
|
|
sscanf(argv[i + 1], "%lf", &minlat);
|
|
if (strcmp(argv[i], "--lonmin") == 0)
|
|
sscanf(argv[i + 1], "%lf", &minlon);
|
|
if (strcmp(argv[i], "--latmax") == 0)
|
|
sscanf(argv[i + 1], "%lf", &maxlat);
|
|
if (strcmp(argv[i], "--lonmax") == 0)
|
|
sscanf(argv[i + 1], "%lf", &maxlon);
|
|
if (strcmp(argv[i], "--out") == 0) {
|
|
filename = argv[i + 1];
|
|
ETOPOMERGING=false;
|
|
}
|
|
if (strcmp(argv[i], "--source") == 0)
|
|
source = argv[i+1];
|
|
}
|
|
|
|
if (sizeof(short int) != 2) {
|
|
cerr << "sizeof(short int) " << sizeof(short int) << endl;
|
|
cerr << "Sorry, can't help you here. Use the source, Luke."
|
|
<< endl;
|
|
exit(1);
|
|
}
|
|
|
|
bool FAILED = false;
|
|
|
|
ifstream etopo2;
|
|
ifstream gtopo30;
|
|
ofstream world3;
|
|
ofstream subset;
|
|
|
|
if (ETOPOMERGING) {
|
|
cout << "Trying to read etopo2.i2 ... ";
|
|
etopo2.open("etopo2.i2", ios::binary);
|
|
if (!etopo2) {
|
|
cout << "failed." << endl;
|
|
FAILED = true;
|
|
} else {
|
|
cout << "ok." << endl;
|
|
}
|
|
}
|
|
|
|
cout << "trying to read " << source << " ... ";
|
|
gtopo30.open(source.c_str(), ios::binary);
|
|
if (!gtopo30) {
|
|
cout << "failed" << endl;
|
|
FAILED = true;
|
|
} else {
|
|
cout << "ok" << endl;
|
|
}
|
|
|
|
if (!ETOPOMERGING) {
|
|
cout << "Trying to open outputfile: " << filename << " ... ";
|
|
subset.open(filename.c_str(), ios::binary);
|
|
if (!subset) {
|
|
cout << "failed." << endl;
|
|
FAILED = true;
|
|
} else {
|
|
cout << "ok." << endl;
|
|
}
|
|
}
|
|
|
|
if (ETOPOMERGING) {
|
|
cout << "Trying to open WORLD3.DEM for writing... ";
|
|
world3.open("WORLD3.DEM", ios::binary);
|
|
if (!world3) {
|
|
cout << "failed." << endl;
|
|
FAILED = true;
|
|
} else {
|
|
cout << "ok." << endl;
|
|
}
|
|
}
|
|
|
|
if (FAILED) failed();
|
|
|
|
int row, col;
|
|
short int alt;
|
|
int i1, i2;
|
|
int offset = 0;
|
|
|
|
// main part of etopo2merging -----------------------------------------------------------------------
|
|
// --------------------------------------------------------------------------------------------------------
|
|
// now read etopo into memory!
|
|
if (ETOPOMERGING) {
|
|
cout << "Reading etopo2 into RAM" << endl;
|
|
for (row = 0; row < 5401; row++) {
|
|
cout << (5401 - row) << " " << '\r' << flush;
|
|
for (col = 0; col < 10801; col++) {
|
|
i1 = etopo2.get();
|
|
i2 = etopo2.get();
|
|
etopo_array[row][col] = 256 * i1 + i2;
|
|
|
|
}
|
|
}
|
|
|
|
cout << "Finished." << endl;
|
|
|
|
cout <<
|
|
"Merging and interpolating... (this may take a little time)" <<
|
|
endl;
|
|
|
|
for (row = 0; row < 21600; row++) {
|
|
cout << (21600 - row) << " " << '\r' << flush;
|
|
for (col = 0; col < 43200; col++) {
|
|
i1 = gtopo30.get();
|
|
i2 = gtopo30.get();
|
|
alt = 256 * i1 + i2;
|
|
if (alt == -9999) { // this is ocean surface, get bathymetry from etopo2
|
|
double lat = 90.0 - (GTOPO30GRID * row);
|
|
double lon = -180.0 + (GTOPO30GRID * col);
|
|
alt = getalt(lat, lon, 0);
|
|
|
|
// Sanity checks
|
|
//if(50473 == offset) {
|
|
// getalt(lat, lon, 1);
|
|
//// exit(1);
|
|
//}
|
|
|
|
// we should not loose the information, that this is covered by water!
|
|
// low areas might have negative elevations!
|
|
// the deepest place on earth is heigher than 16384 (2^14)
|
|
// so we will set the 2^14- bit to indicate water coverage
|
|
|
|
//I havent understand the bit-logic of negative short int..
|
|
|
|
// due to interpolation in the worse resolutioned
|
|
// model we might cat positive values :(
|
|
// they will be forced to 0
|
|
if (alt >= 0) alt =0;
|
|
alt -= 16384;
|
|
}
|
|
|
|
// split alt in big-endian 2 byte and write..
|
|
i2 = (short int) (alt & 255);
|
|
i1 = (short int) (alt & 65280) >> 8;
|
|
|
|
// Sanity checks
|
|
//if(50473 == offset) {
|
|
// printf("i1 = %d, i2 = %d\n", i1, i2);
|
|
//}
|
|
world3.put(i1);
|
|
if (!world3) {
|
|
cerr <<
|
|
"An error occured while writing the output file... (no disc space?)"
|
|
<< endl;
|
|
exit(1);
|
|
}
|
|
|
|
world3.put(i2);
|
|
if (!world3) {
|
|
cerr <<
|
|
"An error occured while writing the output file... (no disc space?)"
|
|
<< endl;
|
|
exit(1);
|
|
}
|
|
offset++;
|
|
}
|
|
}
|
|
}
|
|
|
|
// main part of subset creating .-------------------------------------------------------------------------------
|
|
// -----------------------------------------------------------------------------------------------------------------
|
|
|
|
if (!ETOPOMERGING) {
|
|
cout << "Subset of " << source << " will be created" << endl;
|
|
cout << "Selected limits:" << endl <<
|
|
"Max Lat: " << maxlat << endl <<
|
|
"Min Lat: " << minlat << endl <<
|
|
"Max Lon: " << maxlon << endl <<
|
|
"Min Lon: " << minlon << endl << endl;
|
|
|
|
double grid_lat = 1.0 /120.0;
|
|
double grid_lon = 1.0 / 120.0;
|
|
int downscalefactor = 1;
|
|
|
|
int upsteps = int (((90 - maxlat) / grid_lat)+0.5);
|
|
|
|
int downsteps = int (((90 - minlat) / grid_lat)+0.5);
|
|
|
|
int leftsteps = int (((minlon - (-180)) / grid_lon)+0.5);
|
|
|
|
int rightsteps = int (((maxlon - (-180)) / grid_lon)+0.5);
|
|
|
|
int colums_lon = 43200;
|
|
|
|
double DEM_LAT_MAX = 90.0 -(upsteps * grid_lat);
|
|
double DEM_LAT_MIN = 90.0 - (downsteps * grid_lat);
|
|
double DEM_LON_MAX = -180.0 + (rightsteps * grid_lon);
|
|
double DEM_LON_MIN = -180.0 + (leftsteps * grid_lon);
|
|
|
|
int steps_lat = downsteps - upsteps ; //+ 2 *downscalefactor;
|
|
int steps_lon = rightsteps - leftsteps; // + 2 * downscalefactor;
|
|
|
|
cout << filename << " will be: " << (steps_lat) * (steps_lon )*2 << " bytes." << endl << endl;
|
|
|
|
cout << "DEM_FILE " << filename << endl <<
|
|
"DEM_LAT_MAX " << DEM_LAT_MAX << endl <<
|
|
"DEM_LAT_MIN " << DEM_LAT_MIN << endl <<
|
|
"DEM_LON_MAX " << DEM_LON_MAX << endl <<
|
|
"DEM_LON_MIN " << DEM_LON_MIN << endl <<
|
|
"DEM_ROWS " << steps_lat << endl <<
|
|
"DEM_COLUMNS " << steps_lon << endl;
|
|
|
|
int counter = steps_lat;
|
|
|
|
for (int zeile = 0; zeile < steps_lat; zeile += downscalefactor) {
|
|
cout << counter-- << " " << '\r' << flush;
|
|
for (int spalte = 0; spalte < steps_lon; spalte += downscalefactor) {
|
|
|
|
// fuer 1 byte daten kommentiert ///* */ !?
|
|
|
|
///*
|
|
int verschieber =
|
|
(colums_lon * upsteps + zeile * colums_lon +
|
|
leftsteps + spalte) * 2;
|
|
//*/
|
|
|
|
/*
|
|
int
|
|
verschieber = (colums_lon * upsteps + zeile * colums_lon + leftsteps + spalte);
|
|
*/
|
|
|
|
//cout << "vershieber : " << verschieber << endl;
|
|
gtopo30.seekg(verschieber, ios::beg);
|
|
|
|
///*
|
|
int i1 = gtopo30.get();
|
|
//*/
|
|
|
|
int i2 = gtopo30.get();
|
|
|
|
subset.put(i1);
|
|
subset.put(i2);
|
|
}
|
|
}
|
|
}
|
|
cout << "Finished. Please check the output file." << endl;
|
|
}
|