33#include <gdal_utils.h>
35#include <libdap/DDS.h>
36#include <libdap/ConstraintEvaluator.h>
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>
44#include <dispatch/BESDebug.h>
45#include <dispatch/BESInternalError.h>
46#include <dispatch/TheBESKeys.h>
48#include "FONgTransform.h"
65 d_dest(0), d_dds(dds), d_localfile(localfile),
66 d_geo_transform_set(false), d_width(0.0), d_height(0.0), d_top(0.0), d_left(0.0),
67 d_bottom(0.0), d_right(0.0), d_no_data(0.0), d_no_data_type(none), d_num_bands(0)
69 BESDEBUG(
"fong3",
"dds numvars = " << dds->num_var() << endl);
70 if (localfile.empty())
71 throw BESInternalError(
"Empty local file name passed to constructor", __FILE__, __LINE__);
80 vector<FONgGrid *>::iterator i = d_fong_vars.begin();
81 vector<FONgGrid *>::iterator e = d_fong_vars.end();
91is_convertable_type(
const BaseType *b)
110static FONgGrid *convert(BaseType *v)
114 return new FONgGrid(
static_cast<Grid*
>(v));
117 throw BESInternalError(
"file out GeoTiff, unable to write unknown variable type", __FILE__, __LINE__);
140void FONgTransform::m_scale_data(
double *data)
144 for (
int i = 0; i < width() * height(); ++i)
145 hist.insert(data[i]);
147 BESDEBUG(
"fong3",
"Hist count: " << hist.size() << endl);
149 if (no_data_type() == negative && hist.size() > 1) {
156 set<double>::iterator i = hist.begin();
157 double smallest = *(++i);
158 if (fabs(smallest + no_data()) > 1) {
161 BESDEBUG(
"fong3",
"New no_data value: " << smallest << endl);
163 for (
int i = 0; i < width() * height(); ++i) {
164 if (data[i] <= no_data()) {
171 set<double>::reverse_iterator r = hist.rbegin();
172 double biggest = *(++r);
173 if (fabs(no_data() - biggest) > 1) {
176 BESDEBUG(
"fong3",
"New no_data value: " << biggest << endl);
178 for (
int i = 0; i < width() * height(); ++i) {
179 if (data[i] >= no_data()) {
202 BESDEBUG(
"fong3",
"left: " << d_left <<
", top: " << d_top <<
", right: " << d_right <<
", bottom: " << d_bottom << endl);
203 BESDEBUG(
"fong3",
"width: " << d_width <<
", height: " << d_height << endl);
215 d_gt[1] = (d_right - d_left) / d_width;
216 d_gt[5] = (d_bottom - d_top) / d_height;
218 BESDEBUG(
"fong3",
"gt[0]: " << d_gt[0] <<
", gt[1]: " << d_gt[1] <<
", gt[2]: " << d_gt[2] \
219 <<
", gt[3]: " << d_gt[3] <<
", gt[4]: " << d_gt[4] <<
", gt[5]: " << d_gt[5] << endl);
235bool FONgTransform::effectively_two_D(
FONgGrid *fbtp)
237 if (fbtp->type() == dods_grid_c) {
238 Grid *g =
static_cast<FONgGrid*
>(fbtp)->grid();
240 if (g->get_array()->dimensions() == 2)
244 Array *a = g->get_array();
246 for (Array::Dim_iter d = a->dim_begin(); d != a->dim_end(); ++d) {
247 if (a->dimension_size(d,
true) > 1)
259 if (btp->send_p() && is_convertable_type(btp)) {
260 BESDEBUG(
"fong3",
"converting " << btp->name() << endl);
277 Structure::Vars_iter vi = s->var_begin();
278 while (vi != s->var_end()) {
279 if ((*vi)->send_p() && is_convertable_type(*vi)) {
280 build_delegate(*vi, t);
282 else if ((*vi)->type() == dods_structure_c) {
283 find_vars_helper(
static_cast<Structure*
>(*vi), t);
297 DDS::Vars_iter vi = dds->var_begin();
298 while (vi != dds->var_end()) {
299 BESDEBUG(
"fong3",
"looking at: " << (*vi)->name() <<
" and it is/isn't selected: " << (*vi)->send_p() << endl);
300 if ((*vi)->send_p() && is_convertable_type(*vi)) {
301 BESDEBUG(
"fong3",
"converting " << (*vi)->name() << endl);
302 build_delegate(*vi, t);
304 else if ((*vi)->type() == dods_structure_c) {
305 find_vars_helper(
static_cast<Structure*
>(*vi), t);
322 find_vars(d_dds, *
this);
324 for (
int i = 0; i < num_bands(); ++i)
325 if (!effectively_two_D(var(i)))
326 throw Error(
"GeoTiff responses can consist of two-dimensional variables only; use constraints to reduce the size of Grids as needed.");
328 GDALDriver *Driver = GetGDALDriverManager()->GetDriverByName(
"MEM");
330 throw Error(
"Could not get the MEM driver from/for GDAL: " +
string(CPLGetLastErrorMsg()));
332 char **Metadata = Driver->GetMetadata();
333 if (!CSLFetchBoolean(Metadata, GDAL_DCAP_CREATE, FALSE))
334 throw Error(
"Could not make output format.");
336 BESDEBUG(
"fong3",
"num_bands: " << num_bands() <<
"." << endl);
346 d_dest = Driver->Create(
"in_memory_dataset", width(), height(), num_bands(), GDT_Byte, 0);
348 d_dest = Driver->Create(
"in_memory_dataset", width(), height(), num_bands(), GDT_Float32, 0);
351 throw Error(
"Could not create the geotiff dataset: " +
string(CPLGetLastErrorMsg()));
355 BESDEBUG(
"fong3",
"Made new temp file and set georeferencing (" << num_bands() <<
" vars)." << endl);
357 bool projection_set =
false;
359 for (
int i = 0; i < num_bands(); ++i) {
362 if (!projection_set) {
364 if (d_dest->SetProjection(wkt.c_str()) != CPLE_None)
365 throw Error(
"Could not set the projection: " +
string(CPLGetLastErrorMsg()));
366 projection_set =
true;
371 throw Error(
"In building a multiband response, different bands had different projection information.");
374 GDALRasterBand *band = d_dest->GetRasterBand(i+1);
376 throw Error(
"Could not get the " + long_to_string(i+1) +
"th band: " +
string(CPLGetLastErrorMsg()));
388 vector<double> local_lat;
390 extract_double_array(fbtp->d_lat, local_lat);
395 if (local_lat[0] < local_lat[local_lat.size() - 1]) {
396 BESDEBUG(
"fong3",
"Writing reversed raster. Lat[0] = " << local_lat[0] << endl);
398 for (
int row = 0; row <= height()-1; ++row) {
399 int offsety=height()-row-1;
400 CPLErr error_write = band->RasterIO(GF_Write, 0, offsety, width(), 1, data+(row*width()), width(), 1, GDT_Float64, 0, 0);
401 if (error_write != CPLE_None)
402 throw Error(
"Could not write data for band: " + long_to_string(i + 1) +
": " +
string(CPLGetLastErrorMsg()));
406 BESDEBUG(
"fong3",
"calling band->RasterIO" << endl);
407 CPLErr error = band->RasterIO(GF_Write, 0, 0, width(), height(), data, width(), height(), GDT_Float64, 0, 0);
408 if (error != CPLE_None)
409 throw Error(
"Could not write data for band: " + long_to_string(i+1) +
": " +
string(CPLGetLastErrorMsg()));
422 GDALDataset *tif_dst = 0;
424 Driver = GetGDALDriverManager()->GetDriverByName(
"GTiff");
426 throw Error(
"Could not get driver for GeoTiff: " +
string(CPLGetLastErrorMsg()));
429 char **Metadata = Driver->GetMetadata();
430 if (!CSLFetchBoolean(Metadata, GDAL_DCAP_CREATECOPY, FALSE))
431 BESDEBUG(
"fong",
"Driver does not support dataset creation via 'CreateCopy()'." << endl);
435 char **options = NULL;
436 options = CSLSetNameValue(options,
"PHOTOMETRIC",
"MINISBLACK" );
442 BESDEBUG(
"fong3",
"Before CreateCopy, number of bands: " << d_dest->GetRasterCount() << endl);
446 argv = CSLAddString(argv,
"-scale");
447 GDALTranslateOptions *opts = GDALTranslateOptionsNew(argv, NULL );
448 int usage_error = CE_None;
449 GDALDatasetH dst_handle = GDALTranslate(d_localfile.c_str(), d_dest, opts, &usage_error);
450 if (!dst_handle || usage_error != CE_None) {
451 GDALClose(dst_handle);
452 GDALTranslateOptionsFree(opts);
453 string msg = string(
"Error calling GDAL translate: ") + CPLGetLastErrorMsg();
454 BESDEBUG(
"fong3",
"ERROR transform_to_geotiff(): " << msg << endl);
455 throw BESError(msg, BES_INTERNAL_ERROR, __FILE__, __LINE__);
458 tif_dst = Driver->CreateCopy(d_localfile.c_str(),
reinterpret_cast<GDALDataset*
>(dst_handle), FALSE,
459 options, NULL, NULL);
461 GDALClose(dst_handle);
462 GDALTranslateOptionsFree(opts);
465 throw Error(
"Could not create the GeoTiff dataset: " +
string(CPLGetLastErrorMsg()));
494 find_vars(d_dds, *
this);
496 for (
int i = 0; i < num_bands(); ++i)
497 if (!effectively_two_D(var(i)))
498 throw Error(
"GeoTiff responses can consist of two-dimensional variables only; use constraints to reduce the size of Grids as needed.");
500 GDALDriver *Driver = GetGDALDriverManager()->GetDriverByName(
"MEM");
502 throw Error(
"Could not get the MEM driver from/for GDAL: " +
string(CPLGetLastErrorMsg()));
504 char **Metadata = Driver->GetMetadata();
505 if (!CSLFetchBoolean(Metadata, GDAL_DCAP_CREATE, FALSE))
506 throw Error(
"Driver JP2OpenJPEG does not support dataset creation.");
510 d_dest = Driver->Create(
"in_memory_dataset", width(), height(), num_bands(), GDT_Byte, 0 );
512 throw Error(
"Could not create in-memory dataset: " +
string(CPLGetLastErrorMsg()));
516 BESDEBUG(
"fong3",
"Made new temp file and set georeferencing (" << num_bands() <<
" vars)." << endl);
518 bool projection_set =
false;
520 for (
int i = 0; i < num_bands(); ++i) {
523 if (!projection_set) {
525 if (d_dest->SetProjection(wkt.c_str()) != CPLE_None)
526 throw Error(
"Could not set the projection: " +
string(CPLGetLastErrorMsg()));
527 projection_set =
true;
532 throw Error(
"In building a multiband response, different bands had different projection information.");
535 GDALRasterBand *band = d_dest->GetRasterBand(i+1);
537 throw Error(
"Could not get the " + long_to_string(i+1) +
"th band: " +
string(CPLGetLastErrorMsg()));
546 vector<double> local_lat;
548 extract_double_array(fbtp->d_lat, local_lat);
553 if (local_lat[0] < local_lat[local_lat.size() - 1]) {
554 BESDEBUG(
"fong3",
"Writing reversed raster. Lat[0] = " << local_lat[0] << endl);
556 for (
int row = 0; row <= height()-1; ++row) {
557 int offsety=height()-row-1;
558 CPLErr error_write = band->RasterIO(GF_Write, 0, offsety, width(), 1, data+(row*width()), width(), 1, GDT_Float64, 0, 0);
559 if (error_write != CPLE_None)
560 throw Error(
"Could not write data for band: " + long_to_string(i + 1) +
": " +
string(CPLGetLastErrorMsg()));
564 BESDEBUG(
"fong3",
"calling band->RasterIO" << endl);
565 CPLErr error = band->RasterIO(GF_Write, 0, 0, width(), height(), data, width(), height(), GDT_Float64, 0, 0);
566 if (error != CPLE_None)
567 throw Error(
"Could not write data for band: " + long_to_string(i+1) +
": " +
string(CPLGetLastErrorMsg()));
581 GDALDataset *jpeg_dst = 0;
583 Driver = GetGDALDriverManager()->GetDriverByName(
"JP2OpenJPEG");
585 throw Error(
"Could not get driver for JP2OpenJPEG: " +
string(CPLGetLastErrorMsg()));
588 char **Metadata = Driver->GetMetadata();
589 if (!CSLFetchBoolean(Metadata, GDAL_DCAP_CREATECOPY, FALSE))
590 BESDEBUG(
"fong",
"Driver JP2OpenJPEG does not support dataset creation via 'CreateCopy()'." << endl);
593 char **options = NULL;
594 options = CSLSetNameValue(options,
"CODEC",
"JP2");
595 options = CSLSetNameValue(options,
"GMLJP2",
"YES");
596 options = CSLSetNameValue(options,
"GeoJP2",
"NO");
597 options = CSLSetNameValue(options,
"QUALITY",
"100");
598 options = CSLSetNameValue(options,
"REVERSIBLE",
"YES");
600 BESDEBUG(
"fong3",
"Before JPEG2000 CreateCopy, number of bands: " << d_dest->GetRasterCount() << endl);
604 argv = CSLAddString(argv,
"-scale");
605 GDALTranslateOptions *opts = GDALTranslateOptionsNew(argv, NULL );
606 int usage_error = CE_None;
607 GDALDatasetH dst_h = GDALTranslate(d_localfile.c_str(), d_dest, opts, &usage_error);
608 if (!dst_h || usage_error != CE_None) {
610 GDALTranslateOptionsFree(opts);
611 string msg = string(
"Error calling GDAL translate: ") + CPLGetLastErrorMsg();
612 BESDEBUG(
"fong3",
"ERROR transform_to_jpeg2000(): " << msg << endl);
613 throw BESError(msg, BES_INTERNAL_ERROR, __FILE__, __LINE__);
616 jpeg_dst = Driver->CreateCopy(d_localfile.c_str(),
reinterpret_cast<GDALDataset*
>(dst_h), FALSE,
617 options, NULL, NULL);
620 GDALTranslateOptionsFree(opts);
623 throw Error(
"Could not create the JPEG200 dataset: " +
string(CPLGetLastErrorMsg()));
627 GDALClose (jpeg_dst);
Base exception class for the BES with basic string message.
exception thrown if internal error encountered
A DAP Grid with file out netcdf information included.
string get_projection(libdap::DDS *dds)
Set the projection information For Grids, look for CF information. If it's not present,...
virtual double * get_data()
Get the data values for the band(s). Call must delete.
virtual void extract_coordinates(FONgTransform &t)
Get the GDAL/OGC WKT projection string.
static TheBESKeys * TheKeys()