1 #include <iomanip>
2 #include <sstream>
3
4 #include <simgear/misc/sg_dir.hxx>
5 #include <simgear/debug/logstream.hxx>
6
7 #include <simgear/scene/dem/SGDem.hxx>
8 #include <simgear/scene/dem/SGDemTile.hxx>
9 #include <simgear/scene/dem/SGDemSession.hxx>
10
SGDemTile(const SGPath & levelDir,unsigned wo,unsigned so,int w,int h,int x,int y,int o,bool cache)11 SGDemTile::SGDemTile( const SGPath& levelDir, unsigned wo, unsigned so, int w, int h, int x, int y, int o, bool cache)
12 {
13 // add the tile name to the level path
14 ref_lon = (int)wo/8 - 180;
15 ref_lat = (int)so/8 - 90;
16 path = levelDir / getTileName( ref_lon, ref_lat );
17 width = w;
18 height = h;
19 resx = x;
20 resy = y;
21 overlap = o;
22
23 pixResX = ((double)width/(double)8.0)/(double)(resx-1);
24 pixResY = ((double)height/(double)8.0)/(double)(resy-1);
25
26 if ( cache ) {
27 raster = cacheTile( path );
28 } else {
29 raster = NULL;
30 }
31 }
32
dbgDumpDataset(GDALDataset * poDataset) const33 void SGDemTile::dbgDumpDataset( GDALDataset* poDataset ) const
34 {
35 double adfGeoTransform[6];
36
37 SG_LOG( SG_TERRAIN, SG_INFO, "Driver: " << poDataset->GetDriver()->GetDescription() << "/" << poDataset->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ) );
38 SG_LOG( SG_TERRAIN, SG_INFO, "Size is " << poDataset->GetRasterXSize() << " x " << poDataset->GetRasterYSize() << " x " << poDataset->GetRasterCount() );
39
40 if( poDataset->GetProjectionRef() != NULL ) {
41 SG_LOG( SG_TERRAIN, SG_INFO, "Projection is " << poDataset->GetProjectionRef() );
42 }
43
44 if( poDataset->GetGeoTransform( adfGeoTransform ) == CE_None ) {
45 SG_LOG( SG_TERRAIN, SG_INFO, "Origin = (" << adfGeoTransform[0] << ", " << adfGeoTransform[3] << ")" );
46 SG_LOG( SG_TERRAIN, SG_INFO, "Pixel Size = (" << adfGeoTransform[1] << ", " << adfGeoTransform[5] << ")" );
47 }
48 }
49
dbgDumpBand(GDALRasterBand * poBand) const50 void SGDemTile::dbgDumpBand( GDALRasterBand* poBand ) const
51 {
52 int nBlockXSize, nBlockYSize;
53 int bGotMin, bGotMax;
54 double adfMinMax[2];
55
56 poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
57 SG_LOG( SG_TERRAIN, SG_INFO, "Block=" << nBlockXSize << " x " << nBlockYSize << " Type=" << GDALGetDataTypeName(poBand->GetRasterDataType()) << "ColorInterp=" << GDALGetColorInterpretationName( poBand->GetColorInterpretation()) );
58
59 adfMinMax[0] = poBand->GetMinimum( &bGotMin );
60 adfMinMax[1] = poBand->GetMaximum( &bGotMax );
61 if( ! (bGotMin && bGotMax) ) {
62 GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);
63 }
64 SG_LOG( SG_TERRAIN, SG_INFO, "Min=" << adfMinMax[0] << ", Max=" << adfMinMax[1] );
65 if( poBand->GetOverviewCount() > 0 ) {
66 SG_LOG( SG_TERRAIN, SG_INFO, "Band has " << poBand->GetOverviewCount() << " overviews." );
67 }
68 if( poBand->GetColorTable() != NULL ) {
69 SG_LOG( SG_TERRAIN, SG_INFO, "Band has a color table with " << poBand->GetColorTable()->GetColorEntryCount() << " entries." );
70 }
71 }
72
cacheTile(const SGPath & path)73 unsigned short* SGDemTile::cacheTile( const SGPath& path )
74 {
75 unsigned short* r = NULL;
76
77 GDALDataset* poDataset = NULL;
78 GDALRasterBand* poBand = NULL;
79
80 // check if file exists to supress GDAL errors...
81 if ( path.exists() ) {
82 poDataset = (GDALDataset *)GDALOpen( path.c_str(), GA_ReadOnly );
83 }
84
85 if ( poDataset ) {
86
87 #if DEM_DEBUG
88 dbgDumpDataset( poDataset );
89 #endif
90
91 // read the bands into array
92 unsigned int numRasters = poDataset->GetRasterCount();
93 for ( unsigned int rb=1; rb<=numRasters; rb++ ) {
94 poBand = poDataset->GetRasterBand(rb);
95
96 #if DEM_DEBUG
97 dbgDumpBand( poBand );
98 #endif
99
100 // if this is the raster we're looking for
101 if ( rb == 1 ) {
102 int nXSize = poBand->GetXSize();
103 int nYSize = poBand->GetYSize();
104 // SG_LOG( SG_TERRAIN, SG_INFO, "Reading raster " << rb << " (" << nXSize << "x" << nYSize << ")" );
105
106 r = (unsigned short *)CPLMalloc(sizeof(unsigned short)*nXSize*nYSize);
107 //fprintf( stderr, "reading raster size %d x %d - buffer is %p\n", nXSize, nYSize, r );
108 CPLErr err = poBand->RasterIO( GF_Read, 0, 0, nXSize, nYSize, r, nXSize, nYSize, GDT_UInt16, 0, 0 );
109 if ( err ) {
110 fprintf(stderr, "Error reading raster\n");
111 }
112 }
113 }
114
115 GDALClose( poDataset );
116 }
117
118 return r;
119 }
120
121 // Create a new DEM tile from multiple source tiles ( at an expected lower resolution )
122 // code based on gdalwarp, but with most options removed.
123 // example dgalwarp usage this function is based on:
124 // gdalwarp -r cubic -ts 1201 1201 -te -85.0 32.0 -83.0 34.0 -dstnodata 0 -co "COMPRESS=DEFLATE" temp/N32W085.hgt temp/merged_cubic.tiff
125
SGDemTile(const SGPath & levelDir,unsigned wo,unsigned so,int w,int h,int x,int y,int overlap,const SGDemSession & s,bool & bWritten)126 SGDemTile::SGDemTile( const SGPath& levelDir, unsigned wo, unsigned so, int w, int h, int x, int y, int overlap, const SGDemSession& s, bool& bWritten )
127 {
128 fprintf( stderr, "Writing tile: lon offset is %u, lat offset is %u\n", wo, so );
129
130 // assume failure
131 bWritten = false;
132
133 // add the tile name to the level path
134 ref_lon = (int)wo/8 - 180;
135 ref_lat = (int)so/8 - 90;
136
137 path = levelDir / getTileName( ref_lon, ref_lat );
138 width = w;
139 height = h;
140 resx = x;
141 resy = y;
142
143 raster = NULL;
144
145 // use gdal warp api to generate tile from session
146 char** papszSrcFiles = NULL;
147 char** papszWarpOptions = NULL;
148 char** papszTO = NULL;
149
150 double dfMinX=0.0, dfMinY=0.0, dfMaxX=0.0, dfMaxY=0.0; // -te
151 double overlapw = (double)w/(double)resx * overlap;
152 double overlaph = (double)h/(double)resy * overlap;
153 int nForcePixels = resx+(2*overlap), nForceLines = resy+(2*overlap); // -ts
154
155 // -dstnodata
156 papszWarpOptions = CSLSetNameValue(papszWarpOptions, "INIT_DEST", "0");
157
158 // target extents ( +/- 1 pixel for normals, and skirts )
159 dfMinX = (double)ref_lon - overlapw;
160 dfMinY = (double)ref_lat - overlaph;
161 dfMaxX = (double)ref_lon + w + overlapw;
162 dfMaxY = (double)ref_lat + h + overlaph;
163
164 SG_LOG( SG_TERRAIN, SG_INFO, "overlapw: " << overlapw << " overlapw " << overlapw << " resx " << resx << " resy " << resy << " w " << w << " h " << h );
165 SG_LOG( SG_TERRAIN, SG_INFO, " dfMinX: " << dfMinX << " dfMinY: " << dfMinY << " dfMaxX: " << dfMaxX << " dfMaxY: " << dfMaxY );
166
167 // generate list of source files in session
168 const std::vector<SGDemTileRef>& tiles = s.getTiles();
169 SG_LOG( SG_TERRAIN, SG_INFO, " create tile from " << tiles.size() << " source tiles" );
170
171 for ( unsigned int i=0; i<tiles.size(); i++ ) {
172 // check if the file exists
173 if ( tiles[i]->getPath().exists() ) {
174 papszSrcFiles = CSLAddString( papszSrcFiles, tiles[i]->getPath().c_str() );
175 SG_LOG( SG_TERRAIN, SG_INFO, " Adding tile " << tiles[i]->getPath() );
176 } else {
177 // SG_LOG( SG_TERRAIN, SG_INFO, " tile " << tiles[i]->getPath() << " doesn't exist" );
178 }
179 }
180
181 // create output
182 if ( papszSrcFiles ) {
183 GDALDatasetH hDstDS = createTile( papszSrcFiles, path.c_str(), nForceLines, nForcePixels,
184 dfMinX, dfMinY, dfMaxX, dfMaxY,
185 papszTO );
186
187 if( hDstDS != NULL ) {
188 /* -------------------------------------------------------------------- */
189 /* Loop over all source files, processing each in turn. */
190 /* -------------------------------------------------------------------- */
191 int iSrc;
192 for( iSrc = 0; papszSrcFiles[iSrc] != NULL; iSrc++ )
193 {
194 doWarp( iSrc, papszSrcFiles[iSrc], hDstDS, papszTO, papszWarpOptions );
195 }
196
197 // manually set min/max to all tiles are consistent when viewing in qgis/grass
198 // we need to get the raster band
199 GDALDataset* poDataset = (GDALDataset *)hDstDS;
200 GDALRasterBand* poBand = NULL;
201
202 // read the bands into array
203 unsigned int numRasters = poDataset->GetRasterCount();
204 for ( unsigned int rb=1; rb<=numRasters; rb++ ) {
205 poBand = poDataset->GetRasterBand(rb);
206
207 // if this is the raster we're looking for
208 if ( rb == 1 ) {
209 double pdfMin, pdfMax, pdfMean, pdfStddev;
210
211 poBand->ComputeStatistics(false, &pdfMin, &pdfMax, &pdfMean, &pdfStddev, NULL, NULL );
212
213 fprintf(stderr, "Got band min as %lf, max as %lf\n", pdfMin, pdfMax );
214
215 // force status to sea level / round of mnt everest
216 pdfMin = 0; pdfMax = 9000;
217 poBand->SetStatistics( pdfMin, pdfMax, pdfMean, pdfStddev);
218
219 fprintf(stderr, "Setting band min to %lf, max to %lf\n", pdfMin, pdfMax );
220 }
221 }
222
223 /* -------------------------------------------------------------------- */
224 /* Final Cleanup. */
225 /* -------------------------------------------------------------------- */
226 CPLErrorReset();
227 GDALFlushCache( hDstDS );
228 GDALClose( hDstDS );
229
230 CSLDestroy( papszSrcFiles );
231 CSLDestroy( papszWarpOptions );
232 CSLDestroy( papszTO );
233
234 GDALDumpOpenDatasets( stderr );
235
236 bWritten = true;
237 }
238 }
239 }
240
~SGDemTile()241 SGDemTile::~SGDemTile()
242 {
243 // free the raster
244 if ( raster ) {
245 CPLFree( raster );
246 }
247 }
248
getAlt(const SGGeod & loc) const249 unsigned short SGDemTile::getAlt( const SGGeod& loc ) const
250 {
251 if ( raster ) {
252 bool debug = false;
253
254 // get lon and lat reletive to sw corner
255 double offLon = loc.getLongitudeDeg() - (double)ref_lon;
256 double offLat = loc.getLatitudeDeg() - (double)ref_lat;
257
258 // raster has a 1 pixel border on all sides
259 // take that into account
260 if ( fabs( offLon ) < 0.000001 ) {
261 debug = true;
262 }
263
264 double fractLon = offLon / (double)width;
265 double fractLat = offLat / (double)height;
266
267 int l = (int)( (double)(resy-1) - ((double)(resy-1)*fractLat) + 1 );
268 int p = (int)( (double)(resx-1) * fractLon + 1 );
269
270 if ( debug ) {
271 printf( "SGDemTile::getAlt at %lf,%lf: offLon is %lf, offLat is %lf, width is %d, height is %d, fractLon is %lf, fractLat is %lf, resx is %d, resy is %d, line %d, pixel %d\n",
272 loc.getLongitudeDeg(), loc.getLatitudeDeg(),
273 offLon, offLat, width, height, fractLon, fractLat,
274 resx, resy, l, p );
275 }
276
277 return raster[l*(resx+2)+p];
278 } else {
279 return 0;
280 }
281 }
282
getGeods(unsigned wo,unsigned so,unsigned eo,unsigned no,int grid_width,int grid_height,unsigned subx,unsigned suby,int incw,int inch,::std::vector<SGGeod> & geods,bool Debug1,bool Debug2)283 void SGDemTile::getGeods( unsigned wo, unsigned so, unsigned eo, unsigned no, int grid_width, int grid_height, unsigned subx, unsigned suby, int incw, int inch, ::std::vector<SGGeod>& geods, bool Debug1, bool Debug2 )
284 {
285 // grid width and height include the skirt
286 // sw and ne do not
287 // we need to find the starting and ending l and p;
288 int startl;//, endl;
289 int startp;//, endp;
290 double startlat, startlon;
291 double endlat, endlon;
292
293 startlon = SGDem::offsetToLongitudeDeg(wo) - incw*pixResX;
294 startlat = SGDem::offsetToLatitudeDeg(so) - inch*pixResY;
295
296 endlon = SGDem::offsetToLongitudeDeg(eo) + incw*pixResX;
297 endlat = SGDem::offsetToLatitudeDeg(no) + inch*pixResY;
298
299 // todo : how to calculate 15- from given data
300 if ( Debug1 ) {
301 printf("resx is %d resy is %d incw is %d incy is %d\n", resx, resy, incw, inch );
302 }
303 if ( raster ) {
304 int di, dj;
305 int p, l;
306 double lonPos, latPos;
307 double maxlon = startlon;
308 double maxlat = startlat;
309
310 startl = (resy-1) + overlap - ( suby * (153-3) ) + inch;
311 startp = 0 + overlap + ( subx * (153-3) ) - incw;
312
313 for ( di = 0, p = startp, lonPos = startlon; di < grid_width; di++, p+=incw, lonPos += incw*pixResX ) {
314 maxlon = lonPos;
315
316 for ( dj = 0, l = startl, latPos = startlat; dj < grid_height; dj++, l-=inch, latPos += inch*pixResY ) {
317 maxlat = latPos;
318
319 SGGeod pos = SGGeod::fromDeg( SGMiscd::normalizePeriodic( -180.0, 180.0, lonPos ),
320 SGMiscd::normalizePeriodic( -180.0, 180.0, latPos ) );
321
322 if ( raster ) {
323 pos.setElevationM( raster[l*(resx+(2*overlap))+p] );
324 }
325
326 geods[di*grid_height+dj] = pos;
327 }
328 }
329
330 if ( fabs( endlon - maxlon ) > 0.0001 ) {
331 printf(" tile overlap error %lf : lon is %lf, startlon is %lf, endlon is %lf, maxlon is %lf. grid_width is %d, incw is %d, pixResX is %lf, (grid_width-1)*incw*pixResX is %lf\n",
332 endlon-maxlon, SGDem::offsetToLongitudeDeg(wo), startlon, endlon, maxlon, grid_width, incw, pixResX, (grid_width-1)*incw*pixResX );
333 }
334
335 if ( fabs( endlat - maxlat ) > 0.0001 ) {
336 printf(" tile overlap error %lf : lat is %lf, startlat is %lf, endlat is %lf, maxlat is %lf. grid_height is %d, inch is %d, pixResY is %lf, (grid_height-1)*inch*pixResY is %lf\n",
337 endlat-maxlat, SGDem::offsetToLatitudeDeg(so), startlat, endlat, maxlat, grid_height, inch, pixResY, (grid_height-1)*inch*pixResY );
338 }
339 }
340 }
341
getTileName(int lon,int lat)342 std::string SGDemTile::getTileName( int lon, int lat )
343 {
344 std::stringstream ss;
345
346 if ( lat >= 0 ) {
347 ss << "N" << std::setw(2) << std::setfill('0') << lat;
348 } else {
349 ss << "S" << std::setw(2) << std::setfill('0') << -lat;
350 }
351
352 if ( lon >= 0 ) {
353 ss << "E" << std::setw(3) << std::setfill('0') << lon;
354 } else {
355 ss << "W" << std::setw(3) << std::setfill('0') << -lon;
356 }
357 ss << ".hgt";
358
359 printf("created tile string %s from %d,%d\n", ss.str().c_str(), lon, lat );
360 return ss.str();
361 }
362