bes Updated for version 3.20.10
FONgTransform.cc
1// FONgTransform.cc
2
3// This file is part of BES GDAL File Out Module
4
5// Copyright (c) 2012 OPeNDAP, Inc.
6// Author: James Gallagher <jgallagher@opendap.org>
7//
8// This library is free software; you can redistribute it and/or
9// modify it under the terms of the GNU Lesser General Public
10// License as published by the Free Software Foundation; either
11// version 2.1 of the License, or (at your option) any later version.
12//
13// This library is distributed in the hope that it will be useful,
14// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// Lesser General Public License for more details.
17//
18// You should have received a copy of the GNU Lesser General Public
19// License along with this library; if not, write to the Free Software
20// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21//
22// You can contact University Corporation for Atmospheric Research at
23// 3080 Center Green Drive, Boulder, CO 80301
24
25#include "config.h"
26
27#include <cstdlib>
28
29#include <set>
30
31#include <gdal.h>
32#include <gdal_priv.h>
33#include <gdal_utils.h>
34
35#include <libdap/DDS.h>
36#include <libdap/ConstraintEvaluator.h>
37
38#include <libdap/Structure.h>
39#include <libdap/Array.h>
40#include <libdap/Grid.h>
41#include <libdap/util.h>
42#include <libdap/Error.h>
43
44#include <BESDebug.h>
45#include <BESInternalError.h>
46
47#include "FONgRequestHandler.h"
48#include "FONgTransform.h"
49
50#include "FONgGrid.h"
51
52using namespace std;
53using namespace libdap;
54
65FONgTransform::FONgTransform(DDS *dds, ConstraintEvaluator &/*evaluator*/, const string &localfile) :
66 d_dest(0), d_dds(dds), d_localfile(localfile),
67 d_geo_transform_set(false), d_width(0.0), d_height(0.0), d_top(0.0), d_left(0.0),
68 d_bottom(0.0), d_right(0.0), d_no_data(0.0), d_no_data_type(none), d_num_bands(0)
69{
70 BESDEBUG("fong3", "dds numvars = " << dds->num_var() << endl);
71 if (localfile.empty())
72 throw BESInternalError("Empty local file name passed to constructor", __FILE__, __LINE__);
73}
74
80{
81 vector<FONgGrid *>::iterator i = d_fong_vars.begin();
82 vector<FONgGrid *>::iterator e = d_fong_vars.end();
83 while (i != e) {
84 delete (*i++);
85 }
86}
87
91static bool
92is_convertable_type(const BaseType *b)
93{
94 switch (b->type()) {
95 case dods_grid_c:
96 return true;
97
98 // TODO Add support for Array (?)
99 case dods_array_c:
100 default:
101 return false;
102 }
103}
104
111static FONgGrid *convert(BaseType *v)
112{
113 switch (v->type()) {
114 case dods_grid_c:
115 return new FONgGrid(static_cast<Grid*>(v));
116
117 default:
118 throw BESInternalError("file out GeoTiff, unable to write unknown variable type", __FILE__, __LINE__);
119 }
120}
121
122#if 0
141void FONgTransform::m_scale_data(double *data)
142{
143 // It is an error to call this if no_data_type() is 'none'
144 set<double> hist;
145 for (int i = 0; i < width() * height(); ++i)
146 hist.insert(data[i]);
147
148 BESDEBUG("fong3", "Hist count: " << hist.size() << endl);
149
150 if (no_data_type() == negative && hist.size() > 1) {
151 // Values are sorted by 'weak' <, so this is the smallest value
152 // and ++i is the next smallest value. Assume no_data is the
153 // smallest value in the data set and ++i is the smallest actual
154 // data value. Reset the no_data value to be 1.0 < the smallest
155 // actual value. This makes for a good grayscale photometric
156 // GeoTiff w/o changing the actual data values.
157 set<double>::iterator i = hist.begin();
158 double smallest = *(++i);
159 if (fabs(smallest + no_data()) > 1) {
160 smallest -= 1.0;
161
162 BESDEBUG("fong3", "New no_data value: " << smallest << endl);
163
164 for (int i = 0; i < width() * height(); ++i) {
165 if (data[i] <= no_data()) {
166 data[i] = smallest;
167 }
168 }
169 }
170 }
171 else {
172 set<double>::reverse_iterator r = hist.rbegin();
173 double biggest = *(++r);
174 if (fabs(no_data() - biggest) > 1) {
175 biggest += 1.0;
176
177 BESDEBUG("fong3", "New no_data value: " << biggest << endl);
178
179 for (int i = 0; i < width() * height(); ++i) {
180 if (data[i] >= no_data()) {
181 data[i] = biggest;
182 }
183 }
184 }
185 }
186}
187#endif
188
201{
202
203 BESDEBUG("fong3", "left: " << d_left << ", top: " << d_top << ", right: " << d_right << ", bottom: " << d_bottom << endl);
204 BESDEBUG("fong3", "width: " << d_width << ", height: " << d_height << endl);
205
206 d_gt[0] = d_left; // The leftmost 'x' value, which is longitude
207 d_gt[3] = d_top; // The topmost 'y' value, which is latitude should be max latitude
208
209 // We assume data w/o any rotation (a north-up image)
210 d_gt[2] = 0.0;
211 d_gt[4] = 0.0;
212
213 // Compute the lower left values. Note that wehn GDAL builds the geotiff
214 // output dataset, it correctly inverts the image when the source data has
215 // inverted latitude values.?
216 d_gt[1] = (d_right - d_left) / d_width; // width in pixels; top and bottom in lat
217 d_gt[5] = (d_bottom - d_top) / d_height;
218
219 BESDEBUG("fong3", "gt[0]: " << d_gt[0] << ", gt[1]: " << d_gt[1] << ", gt[2]: " << d_gt[2] \
220 << ", gt[3]: " << d_gt[3] << ", gt[4]: " << d_gt[4] << ", gt[5]: " << d_gt[5] << endl);
221
222 return d_gt;
223}
224
236bool FONgTransform::effectively_two_D(FONgGrid *fbtp)
237{
238 if (fbtp->type() == dods_grid_c) {
239 Grid *g = static_cast<FONgGrid*>(fbtp)->grid();
240
241 if (g->get_array()->dimensions() == 2)
242 return true;
243
244 // count the dimensions with sizes other than 1
245 Array *a = g->get_array();
246 int dims = 0;
247 for (Array::Dim_iter d = a->dim_begin(); d != a->dim_end(); ++d) {
248 if (a->dimension_size(d, true) > 1)
249 ++dims;
250 }
251
252 return dims == 2;
253 }
254
255 return false;
256}
257
258static void build_delegate(BaseType *btp, FONgTransform &t)
259{
260 if (btp->send_p() && is_convertable_type(btp)) {
261 BESDEBUG( "fong3", "converting " << btp->name() << endl);
262
263 // Build the delegate
264 FONgGrid *fb = convert(btp);
265
266 // Get the information needed for the transform.
267 // Note that FONgGrid::extract_coordinates() also pushes the
268 // new FONgBaseType instance onto the FONgTransform's vector of
269 // delagate variable objects.
270 fb->extract_coordinates(t);
271 }
272}
273
274// Helper function to descend into Structures looking for Grids (and Arrays
275// someday).
276static void find_vars_helper(Structure *s, FONgTransform &t)
277{
278 Structure::Vars_iter vi = s->var_begin();
279 while (vi != s->var_end()) {
280 if ((*vi)->send_p() && is_convertable_type(*vi)) {
281 build_delegate(*vi, t);
282 }
283 else if ((*vi)->type() == dods_structure_c) {
284 find_vars_helper(static_cast<Structure*>(*vi), t);
285 }
286
287 ++vi;
288 }
289}
290
291// Helper function to scan the DDS top-level for Grids, ...
292// Note that FONgGrid::extract_coordinates() sets a bunch of
293// values in the FONgBaseType instance _and_ this instance of
294// FONgTransform. One of these is 'num_bands()'. For GeoTiff,
295// num_bands() must be 1. This is tested in transform().
296static void find_vars(DDS *dds, FONgTransform &t)
297{
298 DDS::Vars_iter vi = dds->var_begin();
299 while (vi != dds->var_end()) {
300 BESDEBUG( "fong3", "looking at: " << (*vi)->name() << " and it is/isn't selected: " << (*vi)->send_p() << endl);
301 if ((*vi)->send_p() && is_convertable_type(*vi)) {
302 BESDEBUG( "fong3", "converting " << (*vi)->name() << endl);
303 build_delegate(*vi, t);
304 }
305 else if ((*vi)->type() == dods_structure_c) {
306 find_vars_helper(static_cast<Structure*>(*vi), t);
307 }
308
309 ++vi;
310 }
311}
312
320{
321 // Scan the entire DDS looking for variables that have been 'projected' and
322 // build the delegate objects for them.
323 find_vars(d_dds, *this);
324
325 for (int i = 0; i < num_bands(); ++i)
326 if (!effectively_two_D(var(i)))
327 throw Error("GeoTiff responses can consist of two-dimensional variables only; use constraints to reduce the size of Grids as needed.");
328
329 GDALDriver *Driver = GetGDALDriverManager()->GetDriverByName("MEM");
330 if( Driver == NULL )
331 throw Error("Could not get the MEM driver from/for GDAL: " + string(CPLGetLastErrorMsg()));
332
333 char **Metadata = Driver->GetMetadata();
334 if (!CSLFetchBoolean(Metadata, GDAL_DCAP_CREATE, FALSE))
335 throw Error("Could not make output format.");
336
337 BESDEBUG("fong3", "num_bands: " << num_bands() << "." << endl);
338
339 // Create band in the memory using data type GDT_Byte.
340 // Most image viewers reproduce tiff files with Bits/Sample: 8
341
342 // Make this type depend on the value of a bes.conf parameter.
343 // See FONgRequestHandler.cc and FONgRequestHandler::d_use_byte_for_geotiff_bands.
344 // FIXME This is a hack. But maybe it's good enough?
345 // jhrg 3/20/19
346 if (FONgRequestHandler::get_use_byte_for_geotiff_bands())
347 d_dest = Driver->Create("in_memory_dataset", width(), height(), num_bands(), GDT_Byte, 0/*options*/);
348 else
349 d_dest = Driver->Create("in_memory_dataset", width(), height(), num_bands(), GDT_Float32, 0/*options*/);
350
351 if (!d_dest)
352 throw Error("Could not create the geotiff dataset: " + string(CPLGetLastErrorMsg()));
353
354 d_dest->SetGeoTransform(geo_transform());
355
356 BESDEBUG("fong3", "Made new temp file and set georeferencing (" << num_bands() << " vars)." << endl);
357
358 bool projection_set = false;
359 string wkt = "";
360 for (int i = 0; i < num_bands(); ++i) {
361 FONgGrid *fbtp = var(i);
362
363 if (!projection_set) {
364 wkt = fbtp->get_projection(d_dds);
365 if (d_dest->SetProjection(wkt.c_str()) != CPLE_None)
366 throw Error("Could not set the projection: " + string(CPLGetLastErrorMsg()));
367 projection_set = true;
368 }
369 else {
370 string wkt_i = fbtp->get_projection(d_dds);
371 if (wkt_i != wkt)
372 throw Error("In building a multiband response, different bands had different projection information.");
373 }
374
375 GDALRasterBand *band = d_dest->GetRasterBand(i+1);
376 if (!band)
377 throw Error("Could not get the " + long_to_string(i+1) + "th band: " + string(CPLGetLastErrorMsg()));
378
379 double *data = 0;
380
381 try {
382 // TODO We can read any of the basic DAP2 types and let RasterIO convert it to any other type.
383 // That is, we can read these values in their native type, skipping the conversion here. That
384 // would make this code faster. jhrg 3/20/19
385 data = fbtp->get_data();
386
387 // If the latitude values are inverted, the 0th value will be less than
388 // the last value.
389 vector<double> local_lat;
390 // TODO: Rewrite this to avoid expensive procedure
391 extract_double_array(fbtp->d_lat, local_lat);
392
393 // NB: Here the 'type' value indicates the type of data in the buffer. The
394 // type of the band is set above when the dataset is created.
395
396 if (local_lat[0] < local_lat[local_lat.size() - 1]) {
397 BESDEBUG("fong3", "Writing reversed raster. Lat[0] = " << local_lat[0] << endl);
398 //---------- write reverse raster----------
399 for (int row = 0; row <= height()-1; ++row) {
400 int offsety=height()-row-1;
401 CPLErr error_write = band->RasterIO(GF_Write, 0, offsety, width(), 1, data+(row*width()), width(), 1, GDT_Float64, 0, 0);
402 if (error_write != CPLE_None)
403 throw Error("Could not write data for band: " + long_to_string(i + 1) + ": " + string(CPLGetLastErrorMsg()));
404 }
405 }
406 else {
407 BESDEBUG("fong3", "calling band->RasterIO" << endl);
408 CPLErr error = band->RasterIO(GF_Write, 0, 0, width(), height(), data, width(), height(), GDT_Float64, 0, 0);
409 if (error != CPLE_None)
410 throw Error("Could not write data for band: " + long_to_string(i+1) + ": " + string(CPLGetLastErrorMsg()));
411 }
412
413 delete[] data;
414 }
415 catch (...) {
416 delete[] data;
417 GDALClose(d_dest);
418 throw;
419 }
420 }
421
422 // Now get the GTiff driver and use CreateCopy() on the d_dest "MEM" dataset
423 GDALDataset *tif_dst = 0;
424 try {
425 Driver = GetGDALDriverManager()->GetDriverByName("GTiff");
426 if (Driver == NULL)
427 throw Error("Could not get driver for GeoTiff: " + string(CPLGetLastErrorMsg()));
428
429 // The drivers only support CreateCopy()
430 char **Metadata = Driver->GetMetadata();
431 if (!CSLFetchBoolean(Metadata, GDAL_DCAP_CREATECOPY, FALSE))
432 BESDEBUG("fong", "Driver does not support dataset creation via 'CreateCopy()'." << endl);
433
434 // NB: Changing PHOTOMETIC to MINISWHITE doesn't seem to have any visible affect,
435 // although the resulting files differ. jhrg 11/21/12
436 char **options = NULL;
437 options = CSLSetNameValue(options, "PHOTOMETRIC", "MINISBLACK" ); // The default for GDAL
438 // TODO Adding these options: -co COMPRESS=LZW -co TILED=YES -co INTERLEAVE=BAND
439 // will produce a COG file. The INTERLEAVE=BAND is not strictly needed but would
440 // be nice if there are multiple variables.
441 // See https://geoexamples.com/other/2019/02/08/cog-tutorial.html
442
443 BESDEBUG("fong3", "Before CreateCopy, number of bands: " << d_dest->GetRasterCount() << endl);
444
445 // implementation of gdal_translate -scale to adjust color levels
446 char **argv = NULL;
447 argv = CSLAddString(argv, "-scale");
448 GDALTranslateOptions *opts = GDALTranslateOptionsNew(argv, NULL /*binary options*/);
449 int usage_error = CE_None;
450 GDALDatasetH dst_handle = GDALTranslate(d_localfile.c_str(), d_dest, opts, &usage_error);
451 if (!dst_handle || usage_error != CE_None) {
452 GDALClose(dst_handle);
453 GDALTranslateOptionsFree(opts);
454 string msg = string("Error calling GDAL translate: ") + CPLGetLastErrorMsg();
455 BESDEBUG("fong3", "ERROR transform_to_geotiff(): " << msg << endl);
456 throw BESError(msg, BES_INTERNAL_ERROR, __FILE__, __LINE__);
457 }
458
459 tif_dst = Driver->CreateCopy(d_localfile.c_str(), reinterpret_cast<GDALDataset*>(dst_handle), FALSE/*strict*/,
460 options, NULL/*progress*/, NULL/*progress data*/);
461
462 GDALClose(dst_handle);
463 GDALTranslateOptionsFree(opts);
464
465 if (!tif_dst)
466 throw Error("Could not create the GeoTiff dataset: " + string(CPLGetLastErrorMsg()));
467 }
468 catch (...) {
469 GDALClose(d_dest);
470 GDALClose (tif_dst);
471 throw;
472 }
473 GDALClose(d_dest);
474 GDALClose(tif_dst);
475}
476
492{
493 // Scan the entire DDS looking for variables that have been 'projected' and
494 // build the delegate objects for them.
495 find_vars(d_dds, *this);
496
497 for (int i = 0; i < num_bands(); ++i)
498 if (!effectively_two_D(var(i)))
499 throw Error("GeoTiff responses can consist of two-dimensional variables only; use constraints to reduce the size of Grids as needed.");
500
501 GDALDriver *Driver = GetGDALDriverManager()->GetDriverByName("MEM");
502 if( Driver == NULL )
503 throw Error("Could not get the MEM driver from/for GDAL: " + string(CPLGetLastErrorMsg()));
504
505 char **Metadata = Driver->GetMetadata();
506 if (!CSLFetchBoolean(Metadata, GDAL_DCAP_CREATE, FALSE))
507 throw Error("Driver JP2OpenJPEG does not support dataset creation.");
508
509 // No creation options for a memory dataset
510 // NB: This is where the type of the bands is set. JPEG2000 only supports integer types.
511 d_dest = Driver->Create("in_memory_dataset", width(), height(), num_bands(), GDT_Byte, 0 /*options*/);
512 if (!d_dest)
513 throw Error("Could not create in-memory dataset: " + string(CPLGetLastErrorMsg()));
514
515 d_dest->SetGeoTransform(geo_transform());
516
517 BESDEBUG("fong3", "Made new temp file and set georeferencing (" << num_bands() << " vars)." << endl);
518
519 bool projection_set = false;
520 string wkt = "";
521 for (int i = 0; i < num_bands(); ++i) {
522 FONgGrid *fbtp = var(i);
523
524 if (!projection_set) {
525 wkt = fbtp->get_projection(d_dds);
526 if (d_dest->SetProjection(wkt.c_str()) != CPLE_None)
527 throw Error("Could not set the projection: " + string(CPLGetLastErrorMsg()));
528 projection_set = true;
529 }
530 else {
531 string wkt_i = fbtp->get_projection(d_dds);
532 if (wkt_i != wkt)
533 throw Error("In building a multiband response, different bands had different projection information.");
534 }
535
536 GDALRasterBand *band = d_dest->GetRasterBand(i+1);
537 if (!band)
538 throw Error("Could not get the " + long_to_string(i+1) + "th band: " + string(CPLGetLastErrorMsg()));
539
540 double *data = 0;
541 try {
542 // TODO We can read any of the basic DAP2 types and let RasterIO convert it to any other type.
543 data = fbtp->get_data();
544
545 // If the latitude values are inverted, the 0th value will be less than
546 // the last value.
547 vector<double> local_lat;
548 // TODO: Rewrite this to avoid expensive procedure
549 extract_double_array(fbtp->d_lat, local_lat);
550
551 // NB: Here the 'type' value indicates the type of data in the buffer. The
552 // type of the band is set above when the dataset is created.
553
554 if (local_lat[0] < local_lat[local_lat.size() - 1]) {
555 BESDEBUG("fong3", "Writing reversed raster. Lat[0] = " << local_lat[0] << endl);
556 //---------- write reverse raster----------
557 for (int row = 0; row <= height()-1; ++row) {
558 int offsety=height()-row-1;
559 CPLErr error_write = band->RasterIO(GF_Write, 0, offsety, width(), 1, data+(row*width()), width(), 1, GDT_Float64, 0, 0);
560 if (error_write != CPLE_None)
561 throw Error("Could not write data for band: " + long_to_string(i + 1) + ": " + string(CPLGetLastErrorMsg()));
562 }
563 }
564 else {
565 BESDEBUG("fong3", "calling band->RasterIO" << endl);
566 CPLErr error = band->RasterIO(GF_Write, 0, 0, width(), height(), data, width(), height(), GDT_Float64, 0, 0);
567 if (error != CPLE_None)
568 throw Error("Could not write data for band: " + long_to_string(i+1) + ": " + string(CPLGetLastErrorMsg()));
569 }
570
571 delete[] data;
572
573 }
574 catch (...) {
575 delete[] data;
576 GDALClose(d_dest);
577 throw;
578 }
579 }
580
581 // Now get the OpenJPEG driver and use CreateCopy() on the d_dest "MEM" dataset
582 GDALDataset *jpeg_dst = 0;
583 try {
584 Driver = GetGDALDriverManager()->GetDriverByName("JP2OpenJPEG");
585 if (Driver == NULL)
586 throw Error("Could not get driver for JP2OpenJPEG: " + string(CPLGetLastErrorMsg()));
587
588 // The JPEG2000 drivers only support CreateCopy()
589 char **Metadata = Driver->GetMetadata();
590 if (!CSLFetchBoolean(Metadata, GDAL_DCAP_CREATECOPY, FALSE))
591 BESDEBUG("fong", "Driver JP2OpenJPEG does not support dataset creation via 'CreateCopy()'." << endl);
592 //throw Error("Driver JP2OpenJPEG does not support dataset creation via 'CreateCopy()'.");
593
594 char **options = NULL;
595 options = CSLSetNameValue(options, "CODEC", "JP2");
596 options = CSLSetNameValue(options, "GMLJP2", "YES");
597 options = CSLSetNameValue(options, "GeoJP2", "NO");
598 options = CSLSetNameValue(options, "QUALITY", "100"); // 25 is the default;
599 options = CSLSetNameValue(options, "REVERSIBLE", "YES"); // lossy compression
600
601 BESDEBUG("fong3", "Before JPEG2000 CreateCopy, number of bands: " << d_dest->GetRasterCount() << endl);
602
603 // implementation of gdal_translate -scale to adjust color levels
604 char **argv = NULL;
605 argv = CSLAddString(argv, "-scale");
606 GDALTranslateOptions *opts = GDALTranslateOptionsNew(argv, NULL /*binary options*/);
607 int usage_error = CE_None;
608 GDALDatasetH dst_h = GDALTranslate(d_localfile.c_str(), d_dest, opts, &usage_error);
609 if (!dst_h || usage_error != CE_None) {
610 GDALClose(dst_h);
611 GDALTranslateOptionsFree(opts);
612 string msg = string("Error calling GDAL translate: ") + CPLGetLastErrorMsg();
613 BESDEBUG("fong3", "ERROR transform_to_jpeg2000(): " << msg << endl);
614 throw BESError(msg, BES_INTERNAL_ERROR, __FILE__, __LINE__);
615 }
616
617 jpeg_dst = Driver->CreateCopy(d_localfile.c_str(), reinterpret_cast<GDALDataset*>(dst_h), FALSE/*strict*/,
618 options, NULL/*progress*/, NULL/*progress data*/);
619
620 GDALClose(dst_h);
621 GDALTranslateOptionsFree(opts);
622
623 if (!jpeg_dst)
624 throw Error("Could not create the JPEG200 dataset: " + string(CPLGetLastErrorMsg()));
625 }
626 catch (...) {
627 GDALClose(d_dest);
628 GDALClose (jpeg_dst);
629 throw;
630 }
631
632 GDALClose(d_dest);
633 GDALClose(jpeg_dst);
634}
Abstract exception class for the BES with basic string message.
Definition: BESError.h:58
exception thrown if internal error encountered
A DAP Grid with file out netcdf information included.
Definition: FONgGrid.h:52
string get_projection(libdap::DDS *dds)
Set the projection information For Grids, look for CF information. If it's not present,...
Definition: FONgGrid.cc:296
virtual double * get_data()
Get the data values for the band(s). Call must delete.
Definition: FONgGrid.cc:337
virtual void extract_coordinates(FONgTransform &t)
Get the GDAL/OGC WKT projection string.
Definition: FONgGrid.cc:203
Transformation object that converts an OPeNDAP DataDDS to a GeoTiff file.
Definition: FONgTransform.h:42
FONgTransform(libdap::DDS *dds, libdap::ConstraintEvaluator &evaluator, const string &localfile)
Constructor that creates transformation object from the specified DataDDS object to the specified fil...
virtual ~FONgTransform()
Destructor.
virtual void transform_to_jpeg2000()
Transforms the variables of the DataDDS to a JPEG2000 file.
virtual double * geo_transform()
Build the geotransform array needed by GDAL.
virtual void transform_to_geotiff()
Transforms the variables of the DataDDS to a GeoTiff file.