32 #include <gdal_priv.h>
33 #include <gdal_utils.h>
36 #include <ConstraintEvaluator.h>
38 #include <Structure.h>
45 #include <BESInternalError.h>
47 #include "FONgRequestHandler.h"
48 #include "FONgTransform.h"
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)
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__);
81 vector<FONgGrid *>::iterator i = d_fong_vars.begin();
82 vector<FONgGrid *>::iterator e = d_fong_vars.end();
92 is_convertable_type(
const BaseType *b)
111 static FONgGrid *convert(BaseType *v)
115 return new FONgGrid(
static_cast<Grid*
>(v));
118 throw BESInternalError(
"file out GeoTiff, unable to write unknown variable type", __FILE__, __LINE__);
141 void FONgTransform::m_scale_data(
double *data)
145 for (
int i = 0; i < width() * height(); ++i)
146 hist.insert(data[i]);
148 BESDEBUG(
"fong3",
"Hist count: " << hist.size() << endl);
150 if (no_data_type() == negative && hist.size() > 1) {
157 set<double>::iterator i = hist.begin();
158 double smallest = *(++i);
159 if (fabs(smallest + no_data()) > 1) {
162 BESDEBUG(
"fong3",
"New no_data value: " << smallest << endl);
164 for (
int i = 0; i < width() * height(); ++i) {
165 if (data[i] <= no_data()) {
172 set<double>::reverse_iterator r = hist.rbegin();
173 double biggest = *(++r);
174 if (fabs(no_data() - biggest) > 1) {
177 BESDEBUG(
"fong3",
"New no_data value: " << biggest << endl);
179 for (
int i = 0; i < width() * height(); ++i) {
180 if (data[i] >= no_data()) {
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);
216 d_gt[1] = (d_right - d_left) / d_width;
217 d_gt[5] = (d_bottom - d_top) / d_height;
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);
236 bool FONgTransform::effectively_two_D(
FONgGrid *fbtp)
238 if (fbtp->type() == dods_grid_c) {
239 Grid *g =
static_cast<FONgGrid*
>(fbtp)->grid();
241 if (g->get_array()->dimensions() == 2)
245 Array *a = g->get_array();
247 for (Array::Dim_iter d = a->dim_begin(); d != a->dim_end(); ++d) {
248 if (a->dimension_size(d,
true) > 1)
260 if (btp->send_p() && is_convertable_type(btp)) {
261 BESDEBUG(
"fong3",
"converting " << btp->name() << endl);
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);
283 else if ((*vi)->type() == dods_structure_c) {
284 find_vars_helper(
static_cast<Structure*
>(*vi), t);
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);
305 else if ((*vi)->type() == dods_structure_c) {
306 find_vars_helper(
static_cast<Structure*
>(*vi), t);
323 find_vars(d_dds, *
this);
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.");
329 GDALDriver *Driver = GetGDALDriverManager()->GetDriverByName(
"MEM");
331 throw Error(
"Could not get the MEM driver from/for GDAL: " +
string(CPLGetLastErrorMsg()));
333 char **Metadata = Driver->GetMetadata();
334 if (!CSLFetchBoolean(Metadata, GDAL_DCAP_CREATE, FALSE))
335 throw Error(
"Could not make output format.");
337 BESDEBUG(
"fong3",
"num_bands: " << num_bands() <<
"." << endl);
346 if (FONgRequestHandler::get_use_byte_for_geotiff_bands())
347 d_dest = Driver->Create(
"in_memory_dataset", width(), height(), num_bands(), GDT_Byte, 0);
349 d_dest = Driver->Create(
"in_memory_dataset", width(), height(), num_bands(), GDT_Float32, 0);
352 throw Error(
"Could not create the geotiff dataset: " +
string(CPLGetLastErrorMsg()));
356 BESDEBUG(
"fong3",
"Made new temp file and set georeferencing (" << num_bands() <<
" vars)." << endl);
358 bool projection_set =
false;
360 for (
int i = 0; i < num_bands(); ++i) {
363 if (!projection_set) {
365 if (d_dest->SetProjection(wkt.c_str()) != CPLE_None)
366 throw Error(
"Could not set the projection: " +
string(CPLGetLastErrorMsg()));
367 projection_set =
true;
372 throw Error(
"In building a multiband response, different bands had different projection information.");
375 GDALRasterBand *band = d_dest->GetRasterBand(i+1);
377 throw Error(
"Could not get the " + long_to_string(i+1) +
"th band: " +
string(CPLGetLastErrorMsg()));
389 vector<double> local_lat;
391 extract_double_array(fbtp->d_lat, local_lat);
396 if (local_lat[0] < local_lat[local_lat.size() - 1]) {
397 BESDEBUG(
"fong3",
"Writing reversed raster. Lat[0] = " << local_lat[0] << endl);
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()));
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()));
423 GDALDataset *tif_dst = 0;
425 Driver = GetGDALDriverManager()->GetDriverByName(
"GTiff");
427 throw Error(
"Could not get driver for GeoTiff: " +
string(CPLGetLastErrorMsg()));
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);
436 char **options = NULL;
437 options = CSLSetNameValue(options,
"PHOTOMETRIC",
"MINISBLACK" );
439 BESDEBUG(
"fong3",
"Before CreateCopy, number of bands: " << d_dest->GetRasterCount() << endl);
443 argv = CSLAddString(argv,
"-scale");
444 GDALTranslateOptions *opts = GDALTranslateOptionsNew(argv, NULL );
445 int usage_error = CE_None;
446 GDALDatasetH dst_handle = GDALTranslate(d_localfile.c_str(), d_dest, opts, &usage_error);
447 if (!dst_handle || usage_error != CE_None) {
448 GDALClose(dst_handle);
449 GDALTranslateOptionsFree(opts);
450 string msg = string(
"Error calling GDAL translate: ") + CPLGetLastErrorMsg();
451 BESDEBUG(
"fong3",
"ERROR transform_to_geotiff(): " << msg << endl);
452 throw BESError(msg, BES_INTERNAL_ERROR, __FILE__, __LINE__);
455 tif_dst = Driver->CreateCopy(d_localfile.c_str(),
reinterpret_cast<GDALDataset*
>(dst_handle), FALSE,
456 options, NULL, NULL);
458 GDALClose(dst_handle);
459 GDALTranslateOptionsFree(opts);
462 throw Error(
"Could not create the GeoTiff dataset: " +
string(CPLGetLastErrorMsg()));
491 find_vars(d_dds, *
this);
493 for (
int i = 0; i < num_bands(); ++i)
494 if (!effectively_two_D(var(i)))
495 throw Error(
"GeoTiff responses can consist of two-dimensional variables only; use constraints to reduce the size of Grids as needed.");
497 GDALDriver *Driver = GetGDALDriverManager()->GetDriverByName(
"MEM");
499 throw Error(
"Could not get the MEM driver from/for GDAL: " +
string(CPLGetLastErrorMsg()));
501 char **Metadata = Driver->GetMetadata();
502 if (!CSLFetchBoolean(Metadata, GDAL_DCAP_CREATE, FALSE))
503 throw Error(
"Driver JP2OpenJPEG does not support dataset creation.");
507 d_dest = Driver->Create(
"in_memory_dataset", width(), height(), num_bands(), GDT_Byte, 0 );
509 throw Error(
"Could not create in-memory dataset: " +
string(CPLGetLastErrorMsg()));
513 BESDEBUG(
"fong3",
"Made new temp file and set georeferencing (" << num_bands() <<
" vars)." << endl);
515 bool projection_set =
false;
517 for (
int i = 0; i < num_bands(); ++i) {
520 if (!projection_set) {
522 if (d_dest->SetProjection(wkt.c_str()) != CPLE_None)
523 throw Error(
"Could not set the projection: " +
string(CPLGetLastErrorMsg()));
524 projection_set =
true;
529 throw Error(
"In building a multiband response, different bands had different projection information.");
532 GDALRasterBand *band = d_dest->GetRasterBand(i+1);
534 throw Error(
"Could not get the " + long_to_string(i+1) +
"th band: " +
string(CPLGetLastErrorMsg()));
543 vector<double> local_lat;
545 extract_double_array(fbtp->d_lat, local_lat);
550 if (local_lat[0] < local_lat[local_lat.size() - 1]) {
551 BESDEBUG(
"fong3",
"Writing reversed raster. Lat[0] = " << local_lat[0] << endl);
553 for (
int row = 0; row <= height()-1; ++row) {
554 int offsety=height()-row-1;
555 CPLErr error_write = band->RasterIO(GF_Write, 0, offsety, width(), 1, data+(row*width()), width(), 1, GDT_Float64, 0, 0);
556 if (error_write != CPLE_None)
557 throw Error(
"Could not write data for band: " + long_to_string(i + 1) +
": " +
string(CPLGetLastErrorMsg()));
561 BESDEBUG(
"fong3",
"calling band->RasterIO" << endl);
562 CPLErr error = band->RasterIO(GF_Write, 0, 0, width(), height(), data, width(), height(), GDT_Float64, 0, 0);
563 if (error != CPLE_None)
564 throw Error(
"Could not write data for band: " + long_to_string(i+1) +
": " +
string(CPLGetLastErrorMsg()));
578 GDALDataset *jpeg_dst = 0;
580 Driver = GetGDALDriverManager()->GetDriverByName(
"JP2OpenJPEG");
582 throw Error(
"Could not get driver for JP2OpenJPEG: " +
string(CPLGetLastErrorMsg()));
585 char **Metadata = Driver->GetMetadata();
586 if (!CSLFetchBoolean(Metadata, GDAL_DCAP_CREATECOPY, FALSE))
587 BESDEBUG(
"fong",
"Driver JP2OpenJPEG does not support dataset creation via 'CreateCopy()'." << endl);
590 char **options = NULL;
591 options = CSLSetNameValue(options,
"CODEC",
"JP2");
592 options = CSLSetNameValue(options,
"GMLJP2",
"YES");
593 options = CSLSetNameValue(options,
"GeoJP2",
"NO");
594 options = CSLSetNameValue(options,
"QUALITY",
"100");
595 options = CSLSetNameValue(options,
"REVERSIBLE",
"YES");
597 BESDEBUG(
"fong3",
"Before JPEG2000 CreateCopy, number of bands: " << d_dest->GetRasterCount() << endl);
601 argv = CSLAddString(argv,
"-scale");
602 GDALTranslateOptions *opts = GDALTranslateOptionsNew(argv, NULL );
603 int usage_error = CE_None;
604 GDALDatasetH dst_h = GDALTranslate(d_localfile.c_str(), d_dest, opts, &usage_error);
605 if (!dst_h || usage_error != CE_None) {
607 GDALTranslateOptionsFree(opts);
608 string msg = string(
"Error calling GDAL translate: ") + CPLGetLastErrorMsg();
609 BESDEBUG(
"fong3",
"ERROR transform_to_jpeg2000(): " << msg << endl);
610 throw BESError(msg, BES_INTERNAL_ERROR, __FILE__, __LINE__);
613 jpeg_dst = Driver->CreateCopy(d_localfile.c_str(),
reinterpret_cast<GDALDataset*
>(dst_h), FALSE,
614 options, NULL, NULL);
617 GDALTranslateOptionsFree(opts);
620 throw Error(
"Could not create the JPEG200 dataset: " +
string(CPLGetLastErrorMsg()));
624 GDALClose (jpeg_dst);