bes  Updated for version 3.19.1
scale_util.cc
1 // -*- mode: c++; c-basic-offset:4 -*-
2 
3 // This file is part of libdap, A C++ implementation of the OPeNDAP Data
4 // Access Protocol.
5 
6 // Copyright (c) 2013 OPeNDAP, Inc.
7 // Author: James Gallagher <jgallagher@opendap.org>
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the GNU Lesser General Public
11 // License as published by the Free Software Foundation; either
12 // version 2.1 of the License, or (at your option) any later version.
13 //
14 // This library is distributed in the hope that it will be useful,
15 // but WITHOUT ANY WARRANTY; without even the implied warranty of
16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 // Lesser General Public License for more details.
18 //
19 // You should have received a copy of the GNU Lesser General Public
20 // License along with this library; if not, write to the Free Software
21 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22 //
23 // You can contact OPeNDAP, Inc. at PO Box 112, Saunderstown, RI. 02874-0112.
24 
25 #include "config.h"
26 
27 // #define DODS_DEBUG
28 
29 //#include <float.h>
30 
31 #include <iostream>
32 #include <vector>
33 #include <memory>
34 #include <limits>
35 #include <sstream>
36 #include <cassert>
37 
38 #include <gdal.h>
39 #include <gdal_priv.h>
40 #include <gdal_utils.h>
41 #include <ogr_spatialref.h>
42 #include <gdalwarper.h>
43 
44 #include <Str.h>
45 #include <Float32.h>
46 #include <Int16.h>
47 #include <Array.h>
48 #include <Grid.h>
49 
50 #include <util.h>
51 #include <Error.h>
52 #include <BESDebug.h>
53 #include <BESError.h>
54 #include <BESDapError.h>
55 
56 #include "ScaleGrid.h"
57 
58 #define DEBUG_KEY "geo"
59 
60 using namespace std;
61 using namespace libdap;
62 
63 namespace functions {
64 
65 #if 0
66 
67 Not Used
68 
69 inline static int is_valid_lat(const double lat)
70 {
71  return (lat >= -90 && lat <= 90);
72 }
73 
74 inline static int is_valid_lon(const double lon)
75 {
76  return (lon >= -180 && lon <= 180);
77 }
78 #endif
79 
86 SizeBox get_size_box(Array *x, Array *y)
87 {
88  int src_x_size = x->dimension_size(x->dim_begin(), true);
89  int src_y_size = y->dimension_size(y->dim_begin(), true);
90 
91  return SizeBox(src_x_size, src_y_size);
92 }
93 
101 static bool same_as(const double a, const double b)
102 {
103  // use float's epsilon since double's is too small for these tests
104  return fabs(a - b) <= numeric_limits<float>::epsilon();
105 }
106 
113 bool monotonic_and_uniform(const vector<double> &values, double res)
114 {
115  vector<double>::size_type end_index = values.size() - 1;
116  for (vector<double>::size_type i = 0; i < end_index; ++i) {
117 // BESDEBUG(DEBUG_KEY, "values[" << i+1 << "]: " << values[i+1] << " - values[" << i << "]: " << values[i] << endl);
118 // BESDEBUG(DEBUG_KEY, values[i+1] - values[i] <<" != res: " << res << endl);
119  if (!same_as((values[i+1] - values[i]), res))
120  return false;
121  }
122 
123  return true;
124 }
125 
137 vector<double> get_geotransform_data(Array *x, Array *y, bool test_maps /* default: false*/)
138 {
139 #ifndef NDEBUG
140  test_maps = true;
141 #endif
142 
143  SizeBox size = get_size_box(x, y);
144 
145  y->read();
146  vector<double> y_values(size.y_size);
147  extract_double_array(y, y_values);
148 
149  double res_y = (y_values[y_values.size()-1] - y_values[0]) / (y_values.size() -1);
150 
151  if (test_maps && !monotonic_and_uniform(y_values, res_y)){
152  string msg = "The grids maps/dimensions must be monotonic and uniform (" + y->name() + ").";
153  BESDEBUG(DEBUG_KEY,"ERROR get_geotransform_data(): " << msg << endl);
154  throw BESError(msg,BES_SYNTAX_USER_ERROR,__FILE__,__LINE__);
155  }
156  x->read();
157  vector<double> x_values(size.x_size);
158  extract_double_array(x, x_values);
159 
160  double res_x = (x_values[x_values.size()-1] - x_values[0]) / (x_values.size() -1);
161 
162  if (test_maps && !monotonic_and_uniform(x_values, res_x)){
163  string msg = "The grids maps/dimensions must be monotonic and uniform (" + x->name() + ").";
164  BESDEBUG(DEBUG_KEY,"ERROR get_geotransform_data(): " << msg << endl);
165  throw BESError(msg,BES_SYNTAX_USER_ERROR,__FILE__,__LINE__);
166  }
167  // Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)
168  // Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)
169 
170  // Original comment:
171  // In case of north up images, the GT(2) and GT(4) coefficients are zero,
172  // and the GT(1) is pixel width, and GT(5) is pixel height. The (GT(0),GT(3))
173  // position is the top left corner of the top left pixel of the raster.
174  //
175  // Note that for an inverted dataset, use min_y and res_y for GT(3) and GT(5)
176  //
177  // 10/27/16 I decided to not treat this as lat/lon information but simply
178  // develop a mathematical transform that will be correct for the data as given
179  // _so long as_ the data are monotonic and uniform. jhrg
180 
181  vector<double> geo_transform(6);
182  geo_transform[0] = x_values[0];
183  geo_transform[1] = res_x;
184  geo_transform[2] = 0; // Assumed because the x/y maps are vectors
185  geo_transform[3] = y_values[0];
186  geo_transform[4] = 0;
187  geo_transform[5] = res_y;
188 
189  return geo_transform;
190 }
191 
203 vector<GDAL_GCP> get_gcp_data(Array *x, Array *y, int sample_x, int sample_y)
204 {
205  SizeBox size = get_size_box(x, y);
206 
207  y->read();
208  vector<double> y_values(size.y_size);
209  extract_double_array(y, y_values);
210 
211  x->read();
212  vector<double> x_values(size.x_size);
213  extract_double_array(x, x_values);
214 
215  // Build the GCP list.
216 
217  // Determine how many control points to use. Subset by a factor of M
218  // but never use less than 10% of of the x and y axis values (each)
219  // FIXME Just use given values for now, which will default to 1.
220 
221  // Build the GCP list, sampling as per sample_x and sample_y
222  unsigned long n_gcps = (size.x_size/sample_x) * (size.y_size/sample_y);
223 
224  vector<GDAL_GCP> gcp_list(n_gcps);
225  GDALInitGCPs(n_gcps, &gcp_list[0]); // allocates the 'list'; free with Deinit
226 
227  unsigned long count = 0;
228  for (int i = 0; i < size.x_size; i += sample_x) {
229  // The test for count < n_gcps corrects for discrepancies between integer
230  // division and the simple increment used by the loops. 3/29/17 jhrg
231  for (int j = 0; count < n_gcps && j < size.y_size; j += sample_y) {
232  gcp_list[count].dfGCPLine = j;
233  gcp_list[count].dfGCPPixel = i;
234  gcp_list[count].dfGCPX = x_values[i];
235  gcp_list[count].dfGCPY = y_values[j];
236 
237  ++count;
238  }
239  }
240 
241  return gcp_list;
242 }
243 
244 GDALDataType get_array_type(const Array *a)
245 {
246  switch (const_cast<Array*>(a)->var()->type()) {
247  case dods_byte_c:
248  return GDT_Byte;
249 
250  case dods_uint16_c:
251  return GDT_UInt16;
252 
253  case dods_int16_c:
254  return GDT_Int16;
255 
256  case dods_uint32_c:
257  return GDT_UInt32;
258 
259  case dods_int32_c:
260  return GDT_Int32;
261 
262  case dods_float32_c:
263  return GDT_Float32;
264 
265  case dods_float64_c:
266  return GDT_Float64;
267 
268  // DAP4 support
269  case dods_uint8_c:
270  case dods_int8_c:
271  return GDT_Byte;
272 
273  case dods_uint64_c:
274  case dods_int64_c:
275  default: {
276  string msg = "Cannot perform geo-spatial operations on "
277  + const_cast<Array*>(a)->var()->type_name() + " data.";
278  BESDEBUG(DEBUG_KEY,"ERROR get_array_type(): " << msg << endl);
279  throw BESError(msg,BES_SYNTAX_USER_ERROR,__FILE__,__LINE__);
280 
281  }
282  }
283 }
284 
293 template <typename T>
294 static Array *transfer_values_helper(GDALRasterBand *band, const unsigned long x, const unsigned long y, Array *a)
295 {
296  // get the data
297  vector<T> buf(x * y);
298  CPLErr error = band->RasterIO(GF_Read, 0, 0, x, y, &buf[0], x, y, get_array_type(a), 0, 0);
299 
300  if (error != CPLE_None){
301  string msg = string("Could not extract data for array.") + CPLGetLastErrorMsg();
302  BESDEBUG(DEBUG_KEY,"ERROR transfer_values_helper(): " << msg << endl);
303  throw BESError(msg,BES_SYNTAX_USER_ERROR,__FILE__,__LINE__);
304 }
305  a->set_value(buf, buf.size());
306 
307  return a;
308 }
309 
317 Array *build_array_from_gdal_dataset(GDALDataset *source, const Array *dest)
318 {
319  // Get the GDALDataset size
320  GDALRasterBand *band = source->GetRasterBand(1);
321  unsigned long x = band->GetXSize();
322  unsigned long y = band->GetYSize();
323 
324  // Build a new DAP Array; use the dest Array's element type
325  Array *result = new Array("result", const_cast<Array*>(dest)->var()->ptr_duplicate());
326  result->append_dim(y);
327  result->append_dim(x);
328 
329  // get the data
330  switch (result->var()->type()) {
331  case dods_byte_c:
332  return transfer_values_helper<dods_byte>(source->GetRasterBand(1), x, y, result);
333  break;
334  case dods_uint16_c:
335  return transfer_values_helper<dods_uint16>(source->GetRasterBand(1), x, y, result);
336  break;
337  case dods_int16_c:
338  return transfer_values_helper<dods_int16>(source->GetRasterBand(1), x, y, result);
339  break;
340  case dods_uint32_c:
341  return transfer_values_helper<dods_uint32>(source->GetRasterBand(1), x, y, result);
342  break;
343  case dods_int32_c:
344  return transfer_values_helper<dods_int32>(source->GetRasterBand(1), x, y, result);
345  break;
346  case dods_float32_c:
347  return transfer_values_helper<dods_float32>(source->GetRasterBand(1), x, y, result);
348  break;
349  case dods_float64_c:
350  return transfer_values_helper<dods_float64>(source->GetRasterBand(1), x, y, result);
351  break;
352  case dods_uint8_c:
353  return transfer_values_helper<dods_byte>(source->GetRasterBand(1), x, y, result);
354  break;
355  case dods_int8_c:
356  return transfer_values_helper<dods_int8>(source->GetRasterBand(1), x, y, result);
357  break;
358  case dods_uint64_c:
359  case dods_int64_c:
360  default:
361  string msg = "The source array to a geo-function contained an unsupported numeric type.";
362  BESDEBUG(DEBUG_KEY,"ERROR build_array_from_gdal_dataset(): " << msg << endl);
363  throw BESError(msg,BES_SYNTAX_USER_ERROR,__FILE__,__LINE__);
364 }
365 }
366 
389 void build_maps_from_gdal_dataset(GDALDataset *dst, Array *x_map, Array *y_map, bool name_maps /*default false */)
390 {
391  // get the geo-transform data
392  vector<double> gt(6);
393  dst->GetGeoTransform(&gt[0]);
394 
395  // Get the GDALDataset size
396  GDALRasterBand *band = dst->GetRasterBand(1);
397 
398  // Build Lon map
399  unsigned long x = band->GetXSize(); // x_map_vals
400 
401  if (name_maps) {
402  x_map->append_dim(x, "Longitude");
403  }
404  else {
405  x_map->append_dim(x);
406  }
407 
408  // for each value, use the geo-transform data to compute a value and store it.
409  vector<dods_float32> x_map_vals(x);
410  dods_float32 *cur_x = &x_map_vals[0];
411  dods_float32 *prev_x = cur_x;
412  // x_map_vals[0] = gt[0];
413  *cur_x++ = gt[0];
414  for (unsigned long i = 1; i < x; ++i) {
415  // x_map_vals[i] = gt[0] + i * gt[1];
416  // x_map_vals[i] = x_map_vals[i-1] + gt[1];
417  *cur_x++ = *prev_x++ + gt[1];
418  }
419 
420  x_map->set_value(&x_map_vals[0], x); // copies values to new storage
421 
422  // Build the Lat map
423  unsigned long y = band->GetYSize();
424 
425  if (name_maps) {
426  y_map->append_dim(y, "Latitude");
427  }
428  else {
429  y_map->append_dim(y);
430  }
431 
432  // for each value, use the geo-transform data to compute a value and store it.
433  vector<dods_float32> y_map_vals(y);
434  dods_float32 *cur_y = &y_map_vals[0];
435  dods_float32 *prev_y = cur_y;
436  // y_map_vals[0] = gt[3];
437  *cur_y++ = gt[3];
438  for (unsigned long i = 1; i < y; ++i) {
439  // y_map_vals[i] = gt[3] + i * gt[5];
440  // y_map_vals[i] = y_map_vals[i-1] + gt[5];
441  *cur_y++ = *prev_y++ + gt[5];
442  }
443 
444  y_map->set_value(&y_map_vals[0], y);
445 }
446 
447 void build_maps_from_gdal_dataset_3D(GDALDataset *dst, Array *t_map, Array *x_map, Array *y_map, bool name_maps /*default false */)
448 {
449  // get the geo-transform data
450  vector<double> gt(6);
451  dst->GetGeoTransform(&gt[0]);
452  BESDEBUG(DEBUG_KEY,"build_maps_from_gdal_dataset_3D: Band number: " << dst->GetRasterCount() << endl);
453  //start loop
454  for(int j=1; j<=dst->GetRasterCount(); j++ ){
455  // Get the GDALDataset size
456  GDALRasterBand *band = dst->GetRasterBand(j);
457 
458  // Build Lon map
459  unsigned long x = band->GetXSize(); // x_map_vals
460 
461  if (name_maps) {
462  x_map->append_dim(x, "Longitude");
463  }
464  else {
465  x_map->append_dim(x);
466  }
467 
468  // for each value, use the geo-transform data to compute a value and store it.
469  vector<dods_float32> x_map_vals(x);
470  dods_float32 *cur_x = &x_map_vals[0];
471  dods_float32 *prev_x = cur_x;
472  // x_map_vals[0] = gt[0];
473  *cur_x++ = gt[0];
474  for (unsigned long i = 1; i < x; ++i) {
475  // x_map_vals[i] = gt[0] + i * gt[1];
476  // x_map_vals[i] = x_map_vals[i-1] + gt[1];
477  *cur_x++ = *prev_x++ + gt[1];
478  }
479 
480  x_map->set_value(&x_map_vals[0], x); // copies values to new storage
481 
482  // Build the Lat map
483  unsigned long y = band->GetYSize();
484 
485  if (name_maps) {
486  y_map->append_dim(y, "Latitude");
487  }
488  else {
489  y_map->append_dim(y);
490  }
491 
492  // for each value, use the geo-transform data to compute a value and store it.
493  vector<dods_float32> y_map_vals(y);
494  dods_float32 *cur_y = &y_map_vals[0];
495  dods_float32 *prev_y = cur_y;
496  // y_map_vals[0] = gt[3];
497  *cur_y++ = gt[3];
498  for (unsigned long i = 1; i < y; ++i) {
499  // y_map_vals[i] = gt[3] + i * gt[5];
500  // y_map_vals[i] = y_map_vals[i-1] + gt[5];
501  *cur_y++ = *prev_y++ + gt[5];
502  }
503 
504  y_map->set_value(&y_map_vals[0], y);
505  }//end loop
506 }
507 
516 double get_missing_data_value(const Array *src)
517 {
518  Array *a = const_cast<Array*>(src);
519 
520  // Read this from the 'missing_value' or '_FillValue' attributes
521  string mv_attr = a->get_attr_table().get_attr("missing_value");
522  if (mv_attr.empty()) mv_attr = a->get_attr_table().get_attr("_FillValue");
523 
524  double missing_data = numeric_limits<double>::quiet_NaN();
525  if (!mv_attr.empty()) {
526  char *endptr;
527  missing_data = strtod(mv_attr.c_str(), &endptr);
528  if (missing_data == 0.0 && endptr == mv_attr.c_str())
529  missing_data = numeric_limits<double>::quiet_NaN();
530  }
531 
532  return missing_data;
533 }
534 
541 Array::Dim_iter get_x_dim(const libdap::Array *src)
542 {
543  Array *a = const_cast<Array*>(src);
544  int numDims = a->dimensions();
545  if (numDims < 2) {
546  stringstream ss;
547  ss << "Ouch! Retrieving the 'x' dimension for the array ";
548  a->print_decl(ss, "", false, true, true);
549  ss << " FAILED Because it has less than 2 dimensions" << endl;
550  BESDEBUG(DEBUG_KEY, ss.str());
551  throw BESError(ss.str(),BES_SYNTAX_USER_ERROR,__FILE__,__LINE__);
552  }
553  Array::Dim_iter start = a->dim_begin();
554  Array::Dim_iter xDim = start + numDims - 1;
555 
556  return xDim;
557 }
558 
565 Array::Dim_iter get_y_dim(const libdap::Array *src)
566 {
567  Array *a = const_cast<Array*>(src);
568  int numDims = a->dimensions();
569  if (numDims < 2) {
570  stringstream ss;
571  ss << "Ouch! Retrieving the 'y' dimension for the array ";
572  a->print_decl(ss, "", false, true, true);
573  ss << " FAILED Because it has less than 2 dimensions" << endl;
574  BESDEBUG(DEBUG_KEY, ss.str());
575  throw BESError(ss.str(),BES_SYNTAX_USER_ERROR,__FILE__,__LINE__);
576  }
577  Array::Dim_iter start = a->dim_begin();
578  Array::Dim_iter yDim = start + numDims - 2;
579  return yDim;
580 }
581 
582 
593 bool array_is_effectively_2D(const libdap::Array *src)
594 {
595  Array *a = const_cast<Array*>(src);
596  int numDims = a->dimensions();
597  if (numDims == 2) return true;
598  if (numDims < 2) return false;
599  // numDims more than 2. Last dim should be x
600  Array::Dim_iter xDim = get_x_dim(a);
601  for (Array::Dim_iter thisDim = a->dim_begin(); thisDim < xDim; thisDim++) {
602  unsigned long size = a->dimension_size(thisDim, true);
603  if (size > 1) {
604  return true;
605  }
606  }
607  return false;
608 }
609 
621 void read_band_data(const Array *src, GDALRasterBand* band)
622 {
623  Array *a = const_cast<Array*>(src);
624 
625  if (!array_is_effectively_2D(src)) {
626  stringstream ss;
627  ss << "Cannot perform geo-spatial operations on an Array (";
628  ss << a->name() << ") with " << a->dimensions() << " dimensions.";
629  ss << "Because the constrained shape of the array: ";
630  a->print_decl(ss,"",false,true,true);
631  ss << " is not a two-dimensional array." << endl;
632  BESDEBUG(DEBUG_KEY, ss.str());
633  throw BESError(ss.str(), BES_SYNTAX_USER_ERROR, __FILE__, __LINE__);
634  }
635 
636  // unsigned long x = a->dimension_size(a->dim_begin(), true);
637  // unsigned long y = a->dimension_size(a->dim_begin() + 1, true);
638 
639  unsigned long x = a->dimension_size(get_x_dim(a), true);
640  unsigned long y = a->dimension_size(get_y_dim(a), true);
641 
642  a->read(); // Should this code use intern_data()? jhrg 10/11/16
643 
644  // We may be able to use AddBand() to skip the I/O operation here
645  // For now, we use read() to load the data values and get_buf() to
646  // access a pointer to them.
647  CPLErr error = band->RasterIO(GF_Write, 0, 0, x, y, a->get_buf(), x, y, get_array_type(a), 0, 0);
648 
649  if (error != CPLE_None){
650  string msg = "Could not load data for grid '" + a->name() + "' msg: '" + CPLGetLastErrorMsg() + "'";
651  BESDEBUG(DEBUG_KEY,"ERROR read_band_data(): " << msg << endl);
652  throw BESError(msg,BES_SYNTAX_USER_ERROR,__FILE__,__LINE__);
653  }
654 }
655 
666 void add_band_data(const Array *src, GDALDataset* ds)
667 {
668  Array *a = const_cast<Array*>(src);
669 
670  //assert(a->dimensions() == 2);
671 
672  a->read();
673 
674  // The MEMory driver supports the DATAPOINTER option.
675  char **options = NULL;
676  ostringstream oss;
677  oss << reinterpret_cast<unsigned long>(a->get_buf());
678  options = CSLSetNameValue(options, "DATAPOINTER", oss.str().c_str());
679 
680  CPLErr error = ds->AddBand(get_array_type(a), options);
681 
682  CSLDestroy(options);
683 
684  if (error != CPLE_None){
685  string msg ="Could not add data for grid '" + a->name() + "': " + CPLGetLastErrorMsg();
686  BESDEBUG(DEBUG_KEY,"ERROR add_band_data(): " << msg << endl);
687  throw BESError(msg,BES_INTERNAL_ERROR,__FILE__,__LINE__);
688  }
689 }
690 
691 #define ADD_BAND 0
692 
710 auto_ptr<GDALDataset> build_src_dataset(Array *data, Array *x, Array *y, const string &srs)
711 {
712  GDALDriver *driver = GetGDALDriverManager()->GetDriverByName("MEM");
713  if(!driver){
714  string msg = string("Could not get the Memory driver for GDAL: ") + CPLGetLastErrorMsg();
715  BESDEBUG(DEBUG_KEY, "ERROR build_src_dataset(): " << msg << endl);
716  throw BESError(msg,BES_INTERNAL_ERROR,__FILE__,__LINE__);
717  }
718 
719  SizeBox array_size = get_size_box(x, y);
720 
721  // The MEM driver takes no creation options jhrg 10/6/16
722  auto_ptr<GDALDataset> ds(driver->Create("result", array_size.x_size, array_size.y_size,
723  1 /* nBands*/, get_array_type(data), NULL /* driver_options */));
724 
725 #if ADD_BAND
726  add_band_data(data, ds.get());
727 #endif
728 
729  // Get the one band for this dataset and load it with data
730  GDALRasterBand *band = ds->GetRasterBand(1);
731  if (!band) {
732  string msg = "Could not get the GDAL RasterBand for Array '" + data->name() + "': " + CPLGetLastErrorMsg();
733  BESDEBUG(DEBUG_KEY,"ERROR build_src_dataset(): " << msg << endl);
734  throw BESError(msg,BES_INTERNAL_ERROR,__FILE__,__LINE__);
735  }
736 
737  // Set the no data value here; I'm not sure what the affect of using NaN will be... jhrg 10/11/16
738  double no_data = get_missing_data_value(data);
739  band->SetNoDataValue(no_data);
740 
741 #if !ADD_BAND
742  read_band_data(data, band);
743 #endif
744 
745  vector<double> geo_transform = get_geotransform_data(x, y);
746  ds->SetGeoTransform(&geo_transform[0]);
747 
748  OGRSpatialReference native_srs;
749  if (CE_None != native_srs.SetWellKnownGeogCS(srs.c_str())){
750  string msg = "Could not set '" + srs + "' as the dataset native CRS.";
751  BESDEBUG(DEBUG_KEY,"ERROR build_src_dataset(): " << msg << endl);
752  throw BESError(msg,BES_SYNTAX_USER_ERROR,__FILE__,__LINE__);
753  }
754  // I'm not sure what to do about the Projected Coordinate system. jhrg 10/6/16
755  // native_srs.SetUTM( 11, TRUE );
756 
757  // Connect the SRS/CRS to the GDAL Dataset
758  char *pszSRS_WKT = NULL;
759  native_srs.exportToWkt( &pszSRS_WKT );
760  ds->SetProjection( pszSRS_WKT );
761  CPLFree( pszSRS_WKT );
762 
763  return ds;
764 }
765 
776 auto_ptr<GDALDataset> scale_dataset(auto_ptr<GDALDataset> src, const SizeBox &size, const string &crs /*""*/,
777  const string &interp /*nearest*/)
778 {
779  char **argv = NULL;
780  argv = CSLAddString(argv, "-of"); // output format
781  argv = CSLAddString(argv, "MEM");
782 
783  argv = CSLAddString(argv, "-outsize"); // output size
784  ostringstream oss;
785  oss << size.x_size;
786  argv = CSLAddString(argv, oss.str().c_str()); // size x
787  oss.str("");
788  oss << size.y_size;
789  argv = CSLAddString(argv, oss.str().c_str()); // size y
790 
791  argv = CSLAddString(argv, "-b"); // band number
792  argv = CSLAddString(argv, "1");
793 
794  argv = CSLAddString(argv, "-r"); // resampling
795  argv = CSLAddString(argv, interp.c_str()); // {nearest(default),bilinear,cubic,cubicspline,lanczos,average,mode}
796 
797  if (!crs.empty()) {
798  argv = CSLAddString(argv, "-a_srs"); // dst SRS (WKT or "EPSG:n")
799  argv = CSLAddString(argv, crs.c_str());
800  }
801 
802  if (BESISDEBUG(DEBUG_KEY)) {
803  char **local = argv;
804  while (*local) {
805  BESDEBUG(DEBUG_KEY, "argv: " << *local++ << endl);
806  }
807 
808  }
809 #if 0
810  char **local = argv;
811  while (*local) {
812  cerr << "argv: " << *local++ << endl;
813  }
814 #endif
815 
816  GDALTranslateOptions *options = GDALTranslateOptionsNew(argv, NULL /*binary options*/);
817 
818  int usage_error = CE_None; // result
819  GDALDatasetH dst_handle = GDALTranslate("warped_dst", src.get(), options, &usage_error);
820  if (!dst_handle || usage_error != CE_None) {
821  GDALClose(dst_handle);
822  GDALTranslateOptionsFree(options);
823  string msg = string("Error calling GDAL translate: ") + CPLGetLastErrorMsg();
824  BESDEBUG(DEBUG_KEY, "ERROR scale_dataset(): " << msg << endl);
825  throw BESError(msg, BES_INTERNAL_ERROR, __FILE__, __LINE__);
826  }
827 
828  auto_ptr<GDALDataset> dst(static_cast<GDALDataset*>(dst_handle));
829 
830  GDALTranslateOptionsFree(options);
831 
832  return dst;
833 }
834 
835 
836 auto_ptr<GDALDataset> scale_dataset_3D(auto_ptr<GDALDataset> src, const SizeBox &size, const string &crs /*""*/,
837  const string &interp /*nearest*/)
838 {
839  char **argv = NULL;
840  argv = CSLAddString(argv, "-of"); // output format
841  argv = CSLAddString(argv, "MEM");
842 
843  argv = CSLAddString(argv, "-outsize"); // output size
844  ostringstream oss;
845  oss << size.x_size;
846  argv = CSLAddString(argv, oss.str().c_str()); // size x
847  oss.str("");
848  oss << size.y_size;
849  argv = CSLAddString(argv, oss.str().c_str()); // size y
850 
851  // all bands in src
852  int n_bands = src.get()->GetRasterCount();
853  for(int i=0; i < n_bands; i++){
854  oss.str("");
855  oss << i+1;
856  argv = CSLAddString(argv, "-b");
857  argv = CSLAddString(argv, oss.str().c_str());
858  }
859 
860  argv = CSLAddString(argv, "-r"); // resampling
861  argv = CSLAddString(argv, interp.c_str()); // {nearest(default),bilinear,cubic,cubicspline,lanczos,average,mode}
862 
863  if (!crs.empty()) {
864  argv = CSLAddString(argv, "-a_srs"); // dst SRS (WKT or "EPSG:n")
865  argv = CSLAddString(argv, crs.c_str());
866  }
867 
868  if (BESISDEBUG(DEBUG_KEY)) {
869  char **local = argv;
870  while (*local) {
871  BESDEBUG(DEBUG_KEY, "argv: " << *local++ << endl);
872  }
873 
874  }
875 #if 0
876  char **local = argv;
877  while (*local) {
878  cerr << "argv: " << *local++ << endl;
879  }
880 #endif
881 
882  GDALTranslateOptions *options = GDALTranslateOptionsNew(argv, NULL /*binary options*/);
883  int usage_error = CE_None; // result
884  GDALDatasetH dst_handle = GDALTranslate("warped_dst", src.get(), options, &usage_error);
885  if (!dst_handle || usage_error != CE_None) {
886  GDALClose(dst_handle);
887  GDALTranslateOptionsFree(options);
888  string msg = string("Error calling GDAL translate: ") + CPLGetLastErrorMsg();
889  BESDEBUG(DEBUG_KEY, "ERROR scale_dataset(): " << msg << endl);
890  throw BESError(msg, BES_INTERNAL_ERROR, __FILE__, __LINE__);
891  }
892 
893  auto_ptr<GDALDataset> dst(static_cast<GDALDataset*>(dst_handle));
894 
895  GDALTranslateOptionsFree(options);
896 
897  return dst;
898 }
899 
900 
913 Grid *scale_dap_array(const Array *data, const Array *x, const Array *y, const SizeBox &size,
914  const string &crs, const string &interp)
915 {
916  // Build GDALDataset for Grid g with lon and lat maps as given
917  Array *d = const_cast<Array*>(data);
918 
919  auto_ptr<GDALDataset> src = build_src_dataset(d, const_cast<Array*>(x), const_cast<Array*>(y));
920 
921  // scale to the new size, using optional CRS and interpolation params
922  auto_ptr<GDALDataset> dst = scale_dataset(src, size, crs, interp);
923 
924  // Build a result Grid: extract the data, build the maps and assemble
925  auto_ptr<Array> built_data(build_array_from_gdal_dataset(dst.get(), d));
926 
927  auto_ptr<Array> built_lat(new Array(y->name(), new Float32(y->name())));
928  auto_ptr<Array> built_lon(new Array(x->name(), new Float32(x->name())));
929 
930  build_maps_from_gdal_dataset(dst.get(), built_lon.get(), built_lat.get());
931 
932  auto_ptr<Grid> result(new Grid(d->name()));
933  result->set_array(built_data.release());
934  result->add_map(built_lat.release(), false);
935  result->add_map(built_lon.release(), false);
936 
937  return result.release();
938 }
939 
950 Grid *scale_dap_grid(const Grid *g, const SizeBox &size, const string &crs, const string &interp)
951 {
952  string func(__func__);
953  func += "() - ";
954 
955  if(!g){
956  throw BESError(func+"The Grid object is null.",BES_INTERNAL_ERROR,__FILE__,__LINE__);
957  }
958  // Build GDALDataset for Grid g with lon and lat maps as given
959  Array *data = dynamic_cast<Array*>(const_cast<Grid*>(g)->array_var());
960  if(!data){
961  throw BESError(func+"Unable to obtain data array from Grid '"+g->name()+"'",BES_INTERNAL_ERROR,__FILE__,__LINE__);
962  }
963  // return iteration
964  Grid::Map_riter ritr = const_cast<Grid*>(g)->map_rbegin();
965  Array *x = dynamic_cast<Array*>(*ritr);
966  Array *y = dynamic_cast<Array*>(*++ritr);
967 
968  if(!x || !y){
969  throw BESError(func+"Unable to obtain 2 Map arrays from Grid '"+g->name()+"'",BES_INTERNAL_ERROR,__FILE__,__LINE__);
970  }
971 
972  return scale_dap_array(data, x, y, size, crs, interp);
973 }
974 
975 #define ADD_BAND 0
976 
994 auto_ptr<GDALDataset> build_src_dataset_3D(Array *data, Array *t, Array *x, Array *y, const string &srs)
995 {
996  Array *d = dynamic_cast<Array*>(data);
997 
998  GDALDriver *driver = GetGDALDriverManager()->GetDriverByName("MEM");
999  if(!driver){
1000  string msg = string("Could not get the Memory driver for GDAL: ") + CPLGetLastErrorMsg();
1001  BESDEBUG(DEBUG_KEY, "ERROR build_src_dataset(): " << msg << endl);
1002  throw BESError(msg,BES_INTERNAL_ERROR,__FILE__,__LINE__);
1003  }
1004 
1005  SizeBox array_size = get_size_box(x, y);
1006  int nBands = t->length();
1007  BESDEBUG(DEBUG_KEY, "nBands = " << nBands << endl);
1008  int nBytes = data->prototype()->width();
1009  const int data_size = x->length() * y->length();
1010  unsigned int dsize = data_size * nBytes;
1011 
1012  auto_ptr<GDALDataset> ds(driver->Create("result", array_size.x_size, array_size.y_size, nBands, get_array_type(d),
1013  NULL /* driver_options */));
1014  data->read();
1015  // start band loop
1016  for(int i=1; i<=nBands; i++){
1017 
1018  GDALRasterBand *band = ds->GetRasterBand(i);
1019  if (!band) {
1020  string msg = "Could not get the GDAL RasterBand for Array '" + data->name() + "': " + CPLGetLastErrorMsg();
1021  BESDEBUG(DEBUG_KEY,"ERROR build_src_dataset(): " << msg << endl);
1022  throw BESError(msg,BES_INTERNAL_ERROR,__FILE__,__LINE__);
1023  }
1024 
1025  double no_data = get_missing_data_value(data);
1026  band->SetNoDataValue(no_data);
1027 
1028  #if !ADD_BAND
1029  CPLErr error = band->RasterIO(GF_Write, 0, 0, x->length(), y->length(), data->get_buf() + dsize*(i-1), x->length(), y->length(), get_array_type(data), 0, 0);
1030  if (error != CPLE_None)
1031  throw Error("Could not write data for band: " + long_to_string(i) + ": " + string(CPLGetLastErrorMsg()));
1032  #endif
1033 
1034 
1035  } // end band loop
1036  vector<double> geo_transform = get_geotransform_data(x, y);
1037  ds->SetGeoTransform(&geo_transform[0]);
1038 
1039  OGRSpatialReference native_srs;
1040  if (CE_None != native_srs.SetWellKnownGeogCS(srs.c_str())){
1041  string msg = "Could not set '" + srs + "' as the dataset native CRS.";
1042  BESDEBUG(DEBUG_KEY,"ERROR build_src_dataset(): " << msg << endl);
1043  throw BESError(msg,BES_SYNTAX_USER_ERROR,__FILE__,__LINE__);
1044  }
1045 
1046  // Connect the SRS/CRS to the GDAL Dataset
1047  char *pszSRS_WKT = NULL;
1048  native_srs.exportToWkt( &pszSRS_WKT );
1049  ds->SetProjection( pszSRS_WKT );
1050  CPLFree( pszSRS_WKT );
1051 
1052  return ds;
1053 }
1054 
1068 Grid *scale_dap_array_3D(const Array *data, const Array *t, const Array *x, const Array *y, const SizeBox &size,
1069  const string &crs, const string &interp)
1070 {
1071  // Build GDALDataset for Grid g with lon and lat maps as given
1072  Array *d = const_cast<Array*>(data);
1073 
1074  auto_ptr<GDALDataset> src = build_src_dataset_3D(d, const_cast<Array*>(t), const_cast<Array*>(x), const_cast<Array*>(y));
1075 
1076  // scale to the new size, using optional CRS and interpolation params
1077 
1078  auto_ptr<GDALDataset> dst = scale_dataset_3D(src, size, crs, interp);
1079 
1080  // Build a result Grid: extract the data, build the maps and assemble
1081  auto_ptr<Array> built_data(build_array_from_gdal_dataset(dst.get(), d));
1082  auto_ptr<Array> built_time(new Array(t->name(), new Float32(t->name())));
1083  auto_ptr<Array> built_lat(new Array(y->name(), new Float32(y->name())));
1084  auto_ptr<Array> built_lon(new Array(x->name(), new Float32(x->name())));
1085 
1086  build_maps_from_gdal_dataset_3D(dst.get(), built_time.get(), built_lon.get(), built_lat.get());
1087 
1088  auto_ptr<Grid> result(new Grid(d->name()));
1089  result->set_array(built_data.release());
1090  result->add_map(built_time.release(), false);
1091  result->add_map(built_lat.release(), false);
1092  result->add_map(built_lon.release(), false);
1093  BESDEBUG(DEBUG_KEY,"result length: " << result.release()->get_array()->dimensions() << endl);
1094  return result.release();
1095 }
1096 
1097 }
1098  // namespace libdap
STL namespace.
Abstract exception class for the BES with basic string message.
Definition: BESError.h:56