bes  Updated for version 3.20.6
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 
310 #if 0
311 
319 template <typename T>
320 static Vector *transfer_values_helper_vec(GDALRasterBand *band, const int x, const int y, Vector *a)
321 {
322  // get the data
323  vector<T> buf(x * y);
324  CPLErr error = band->RasterIO(GF_Read, 0, 0, x, y, &buf[0], x, y, a->type(), 0, 0);
325 
326  if (error != CPLE_None){
327  string msg = string("Could not extract data for array.") + CPLGetLastErrorMsg();
328  BESDEBUG(DEBUG_KEY,"ERROR transfer_values_helper(): " << msg << endl);
329  throw BESError(msg,BES_SYNTAX_USER_ERROR,__FILE__,__LINE__);
330  }
331  a->set_value(buf, buf.size());
332  return a;
333 }
334 #endif
335 
343 Array *build_array_from_gdal_dataset(GDALDataset *source, const Array *dest)
344 {
345  // Get the GDALDataset size
346  GDALRasterBand *band = source->GetRasterBand(1);
347  unsigned long x = band->GetXSize();
348  unsigned long y = band->GetYSize();
349 
350  // Build a new DAP Array; use the dest Array's element type
351  Array *result = new Array("result", const_cast<Array*>(dest)->var()->ptr_duplicate());
352  result->append_dim(y);
353  result->append_dim(x);
354 
355  // get the data
356  switch (result->var()->type()) {
357  case dods_byte_c:
358  return transfer_values_helper<dods_byte>(source->GetRasterBand(1), x, y, result);
359  break;
360  case dods_uint16_c:
361  return transfer_values_helper<dods_uint16>(source->GetRasterBand(1), x, y, result);
362  break;
363  case dods_int16_c:
364  return transfer_values_helper<dods_int16>(source->GetRasterBand(1), x, y, result);
365  break;
366  case dods_uint32_c:
367  return transfer_values_helper<dods_uint32>(source->GetRasterBand(1), x, y, result);
368  break;
369  case dods_int32_c:
370  return transfer_values_helper<dods_int32>(source->GetRasterBand(1), x, y, result);
371  break;
372  case dods_float32_c:
373  return transfer_values_helper<dods_float32>(source->GetRasterBand(1), x, y, result);
374  break;
375  case dods_float64_c:
376  return transfer_values_helper<dods_float64>(source->GetRasterBand(1), x, y, result);
377  break;
378  case dods_uint8_c:
379  return transfer_values_helper<dods_byte>(source->GetRasterBand(1), x, y, result);
380  break;
381  case dods_int8_c:
382  return transfer_values_helper<dods_int8>(source->GetRasterBand(1), x, y, result);
383  break;
384  case dods_uint64_c:
385  case dods_int64_c:
386  default:
387  string msg = "The source array to a geo-function contained an unsupported numeric type.";
388  BESDEBUG(DEBUG_KEY,"ERROR build_array_from_gdal_dataset(): " << msg << endl);
389  throw BESError(msg,BES_SYNTAX_USER_ERROR,__FILE__,__LINE__);
390 }
391 }
392 
400 Array *build_array_from_gdal_dataset_3D(GDALDataset *source3D, const Array *dest){
401 
402  // DAP array result
403  int t_size = source3D->GetRasterCount();
404  int x_size = source3D->GetRasterBand(1)->GetXSize();
405  int y_size = source3D->GetRasterBand(1)->GetYSize();
406  Array *result = new Array("result", const_cast<Array*>(dest)->var()->ptr_duplicate());
407  result->append_dim(t_size);
408  result->append_dim(y_size);
409  result->append_dim(x_size);
410 
411  vector<dods_float32> res(t_size * x_size * y_size);
412  for(int i=0; i< t_size; i++) {
413  // get the data
414  GDALRasterBand *band = source3D->GetRasterBand(i+1);
415  if (!band)
416  throw Error(string("Could not get the GDALRasterBand for the GDALDataset: ") + CPLGetLastErrorMsg());
417  vector<double> gt(6);
418  source3D->GetGeoTransform(&gt[0]);
419  // Extract data from band
420  vector<dods_float32> buf(x_size * y_size);
421  CPLErr error = band->RasterIO(GF_Read, 0, 0, x_size, y_size, &buf[0], x_size, y_size, get_array_type(dest),
422  0, 0);
423  if (error != CPLE_None)
424  throw Error(string("Could not extract data for translated GDAL Dataset.") + CPLGetLastErrorMsg());
425  if(i==0){
426  res =buf;
427  }else{
428  res.insert(res.end(), buf.begin(), buf.end());
429  }
430  }
431  result->set_value(res, res.size());
432  return result;
433 }
434 
435 
458 void build_maps_from_gdal_dataset(GDALDataset *dst, Array *x_map, Array *y_map, bool name_maps /*default false */)
459 {
460  // get the geo-transform data
461  vector<double> gt(6);
462  dst->GetGeoTransform(&gt[0]);
463 
464  // Get the GDALDataset size
465  GDALRasterBand *band = dst->GetRasterBand(1);
466 
467  // Build Lon map
468  unsigned long x = band->GetXSize(); // x_map_vals
469 
470  if (name_maps) {
471  x_map->append_dim(x, "Longitude");
472  }
473  else {
474  x_map->append_dim(x);
475  }
476 
477  // for each value, use the geo-transform data to compute a value and store it.
478  vector<dods_float32> x_map_vals(x);
479  dods_float32 *cur_x = &x_map_vals[0];
480  dods_float32 *prev_x = cur_x;
481  // x_map_vals[0] = gt[0];
482  *cur_x++ = gt[0];
483  for (unsigned long i = 1; i < x; ++i) {
484  // x_map_vals[i] = gt[0] + i * gt[1];
485  // x_map_vals[i] = x_map_vals[i-1] + gt[1];
486  *cur_x++ = *prev_x++ + gt[1];
487  }
488 
489  x_map->set_value(&x_map_vals[0], x); // copies values to new storage
490 
491  // Build the Lat map
492  unsigned long y = band->GetYSize();
493 
494  if (name_maps) {
495  y_map->append_dim(y, "Latitude");
496  }
497  else {
498  y_map->append_dim(y);
499  }
500 
501  // for each value, use the geo-transform data to compute a value and store it.
502  vector<dods_float32> y_map_vals(y);
503  dods_float32 *cur_y = &y_map_vals[0];
504  dods_float32 *prev_y = cur_y;
505  // y_map_vals[0] = gt[3];
506  *cur_y++ = gt[3];
507  for (unsigned long i = 1; i < y; ++i) {
508  // y_map_vals[i] = gt[3] + i * gt[5];
509  // y_map_vals[i] = y_map_vals[i-1] + gt[5];
510  *cur_y++ = *prev_y++ + gt[5];
511  }
512 
513  y_map->set_value(&y_map_vals[0], y);
514 }
515 
539 void build_maps_from_gdal_dataset_3D(GDALDataset *dst, Array *t, Array *t_map, Array *x_map, Array *y_map, bool name_maps /*default false */)
540 {
541  // get the geo-transform data
542  vector<double> gt(6);
543  dst->GetGeoTransform(&gt[0]);
544 
545  // Get the GDALDataset size
546  GDALRasterBand *band = dst->GetRasterBand(1);
547 
548  // Build Time map
549  int t_size = t->length();
550 
551  if (name_maps) {
552  t_map->append_dim(t_size, "Time");
553  }
554  else {
555  t_map->append_dim(t_size);
556  }
557 
558  vector<dods_float32> t_buf(t_size);
559  t->value(&t_buf[0]);
560 
561  t_map->set_value(&t_buf[0], t_size);
562 
563  // Build Lon map
564  unsigned long x = band->GetXSize(); // x_map_vals
565 
566  if (name_maps) {
567  x_map->append_dim(x, "Longitude");
568  }
569  else {
570  x_map->append_dim(x);
571  }
572 
573  // for each value, use the geo-transform data to compute a value and store it.
574  vector<dods_float32> x_map_vals(x);
575  dods_float32 *cur_x = &x_map_vals[0];
576  dods_float32 *prev_x = cur_x;
577  // x_map_vals[0] = gt[0];
578  *cur_x++ = gt[0];
579  for (unsigned long i = 1; i < x; ++i) {
580  // x_map_vals[i] = gt[0] + i * gt[1];
581  // x_map_vals[i] = x_map_vals[i-1] + gt[1];
582  *cur_x++ = *prev_x++ + gt[1];
583  }
584 
585  x_map->set_value(&x_map_vals[0], x); // copies values to new storage
586 
587  // Build the Lat map
588  unsigned long y = band->GetYSize();
589 
590  if (name_maps) {
591  y_map->append_dim(y, "Latitude");
592  }
593  else {
594  y_map->append_dim(y);
595  }
596 
597  // for each value, use the geo-transform data to compute a value and store it.
598  vector<dods_float32> y_map_vals(y);
599  dods_float32 *cur_y = &y_map_vals[0];
600  dods_float32 *prev_y = cur_y;
601  // y_map_vals[0] = gt[3];
602  *cur_y++ = gt[3];
603  for (unsigned long i = 1; i < y; ++i) {
604  // y_map_vals[i] = gt[3] + i * gt[5];
605  // y_map_vals[i] = y_map_vals[i-1] + gt[5];
606  *cur_y++ = *prev_y++ + gt[5];
607  }
608 
609  y_map->set_value(&y_map_vals[0], y);
610 }
611 
620 double get_missing_data_value(const Array *src)
621 {
622  Array *a = const_cast<Array*>(src);
623 
624  // Read this from the 'missing_value' or '_FillValue' attributes
625  string mv_attr = a->get_attr_table().get_attr("missing_value");
626  if (mv_attr.empty()) mv_attr = a->get_attr_table().get_attr("_FillValue");
627 
628  double missing_data = numeric_limits<double>::quiet_NaN();
629  if (!mv_attr.empty()) {
630  char *endptr;
631  missing_data = strtod(mv_attr.c_str(), &endptr);
632  if (missing_data == 0.0 && endptr == mv_attr.c_str())
633  missing_data = numeric_limits<double>::quiet_NaN();
634  }
635 
636  return missing_data;
637 }
638 
645 Array::Dim_iter get_x_dim(const libdap::Array *src)
646 {
647  Array *a = const_cast<Array*>(src);
648  int numDims = a->dimensions();
649  if (numDims < 2) {
650  stringstream ss;
651  ss << "Ouch! Retrieving the 'x' dimension for the array ";
652  a->print_decl(ss, "", false, true, true);
653  ss << " FAILED Because it has less than 2 dimensions" << endl;
654  BESDEBUG(DEBUG_KEY, ss.str());
655  throw BESError(ss.str(),BES_SYNTAX_USER_ERROR,__FILE__,__LINE__);
656  }
657  Array::Dim_iter start = a->dim_begin();
658  Array::Dim_iter xDim = start + numDims - 1;
659 
660  return xDim;
661 }
662 
669 Array::Dim_iter get_y_dim(const libdap::Array *src)
670 {
671  Array *a = const_cast<Array*>(src);
672  int numDims = a->dimensions();
673  if (numDims < 2) {
674  stringstream ss;
675  ss << "Ouch! Retrieving the 'y' dimension for the array ";
676  a->print_decl(ss, "", false, true, true);
677  ss << " FAILED Because it has less than 2 dimensions" << endl;
678  BESDEBUG(DEBUG_KEY, ss.str());
679  throw BESError(ss.str(),BES_SYNTAX_USER_ERROR,__FILE__,__LINE__);
680  }
681  Array::Dim_iter start = a->dim_begin();
682  Array::Dim_iter yDim = start + numDims - 2;
683  return yDim;
684 }
685 
686 
697 bool array_is_effectively_2D(const libdap::Array *src)
698 {
699  Array *a = const_cast<Array*>(src);
700  int numDims = a->dimensions();
701  if (numDims == 2) return true;
702  if (numDims < 2) return false;
703  // numDims more than 2. Last dim should be x
704  Array::Dim_iter xDim = get_x_dim(a);
705  for (Array::Dim_iter thisDim = a->dim_begin(); thisDim < xDim; thisDim++) {
706  unsigned long size = a->dimension_size(thisDim, true);
707  if (size > 1) {
708  return true;
709  }
710  }
711  return false;
712 }
713 
725 void read_band_data(const Array *src, GDALRasterBand* band)
726 {
727  Array *a = const_cast<Array*>(src);
728 
729  if (!array_is_effectively_2D(src)) {
730  stringstream ss;
731  ss << "Cannot perform geo-spatial operations on an Array (";
732  ss << a->name() << ") with " << a->dimensions() << " dimensions.";
733  ss << "Because the constrained shape of the array: ";
734  a->print_decl(ss,"",false,true,true);
735  ss << " is not a two-dimensional array." << endl;
736  BESDEBUG(DEBUG_KEY, ss.str());
737  throw BESError(ss.str(), BES_SYNTAX_USER_ERROR, __FILE__, __LINE__);
738  }
739 
740  // unsigned long x = a->dimension_size(a->dim_begin(), true);
741  // unsigned long y = a->dimension_size(a->dim_begin() + 1, true);
742 
743  unsigned long x = a->dimension_size(get_x_dim(a), true);
744  unsigned long y = a->dimension_size(get_y_dim(a), true);
745 
746  a->read(); // Should this code use intern_data()? jhrg 10/11/16
747 
748  // We may be able to use AddBand() to skip the I/O operation here
749  // For now, we use read() to load the data values and get_buf() to
750  // access a pointer to them.
751  CPLErr error = band->RasterIO(GF_Write, 0, 0, x, y, a->get_buf(), x, y, get_array_type(a), 0, 0);
752 
753  if (error != CPLE_None){
754  string msg = "Could not load data for grid '" + a->name() + "' msg: '" + CPLGetLastErrorMsg() + "'";
755  BESDEBUG(DEBUG_KEY,"ERROR read_band_data(): " << msg << endl);
756  throw BESError(msg,BES_SYNTAX_USER_ERROR,__FILE__,__LINE__);
757  }
758 }
759 
770 void add_band_data(const Array *src, GDALDataset* ds)
771 {
772  Array *a = const_cast<Array*>(src);
773 
774  //assert(a->dimensions() == 2);
775 
776  a->read();
777 
778  // The MEMory driver supports the DATAPOINTER option.
779  char **options = NULL;
780  ostringstream oss;
781  oss << reinterpret_cast<unsigned long>(a->get_buf());
782  options = CSLSetNameValue(options, "DATAPOINTER", oss.str().c_str());
783 
784  CPLErr error = ds->AddBand(get_array_type(a), options);
785 
786  CSLDestroy(options);
787 
788  if (error != CPLE_None){
789  string msg ="Could not add data for grid '" + a->name() + "': " + CPLGetLastErrorMsg();
790  BESDEBUG(DEBUG_KEY,"ERROR add_band_data(): " << msg << endl);
791  throw BESError(msg,BES_INTERNAL_ERROR,__FILE__,__LINE__);
792  }
793 }
794 
795 #define ADD_BAND 0
796 
814 unique_ptr<GDALDataset> build_src_dataset(Array *data, Array *x, Array *y, const string &srs)
815 {
816  GDALDriver *driver = GetGDALDriverManager()->GetDriverByName("MEM");
817  if(!driver){
818  string msg = string("Could not get the Memory driver for GDAL: ") + CPLGetLastErrorMsg();
819  BESDEBUG(DEBUG_KEY, "ERROR build_src_dataset(): " << msg << endl);
820  throw BESError(msg,BES_INTERNAL_ERROR,__FILE__,__LINE__);
821  }
822 
823  SizeBox array_size = get_size_box(x, y);
824 
825  // The MEM driver takes no creation options jhrg 10/6/16
826  unique_ptr<GDALDataset> ds(driver->Create("result", array_size.x_size, array_size.y_size,
827  1 /* nBands*/, get_array_type(data), NULL /* driver_options */));
828 
829 #if ADD_BAND
830  add_band_data(data, ds.get());
831 #endif
832 
833  // Get the one band for this dataset and load it with data
834  GDALRasterBand *band = ds->GetRasterBand(1);
835  if (!band) {
836  string msg = "Could not get the GDAL RasterBand for Array '" + data->name() + "': " + CPLGetLastErrorMsg();
837  BESDEBUG(DEBUG_KEY,"ERROR build_src_dataset(): " << msg << endl);
838  throw BESError(msg,BES_INTERNAL_ERROR,__FILE__,__LINE__);
839  }
840 
841  // Set the no data value here; I'm not sure what the affect of using NaN will be... jhrg 10/11/16
842  double no_data = get_missing_data_value(data);
843  band->SetNoDataValue(no_data);
844 
845 #if !ADD_BAND
846  read_band_data(data, band);
847 #endif
848 
849  vector<double> geo_transform = get_geotransform_data(x, y);
850  ds->SetGeoTransform(&geo_transform[0]);
851 
852  OGRSpatialReference native_srs;
853  if (CE_None != native_srs.SetWellKnownGeogCS(srs.c_str())){
854  string msg = "Could not set '" + srs + "' as the dataset native CRS.";
855  BESDEBUG(DEBUG_KEY,"ERROR build_src_dataset(): " << msg << endl);
856  throw BESError(msg,BES_SYNTAX_USER_ERROR,__FILE__,__LINE__);
857  }
858  // I'm not sure what to do about the Projected Coordinate system. jhrg 10/6/16
859  // native_srs.SetUTM( 11, TRUE );
860 
861  // Connect the SRS/CRS to the GDAL Dataset
862  char *pszSRS_WKT = NULL;
863  native_srs.exportToWkt( &pszSRS_WKT );
864  ds->SetProjection( pszSRS_WKT );
865  CPLFree( pszSRS_WKT );
866 
867  return ds;
868 }
869 
880 unique_ptr<GDALDataset> scale_dataset(unique_ptr<GDALDataset>& src, const SizeBox &size, const string &crs /*""*/,
881  const string &interp /*nearest*/)
882 {
883  char **argv = NULL;
884  argv = CSLAddString(argv, "-of"); // output format
885  argv = CSLAddString(argv, "MEM");
886 
887  argv = CSLAddString(argv, "-outsize"); // output size
888  ostringstream oss;
889  oss << size.x_size;
890  argv = CSLAddString(argv, oss.str().c_str()); // size x
891  oss.str("");
892  oss << size.y_size;
893  argv = CSLAddString(argv, oss.str().c_str()); // size y
894 
895  argv = CSLAddString(argv, "-b"); // band number
896  argv = CSLAddString(argv, "1");
897 
898  argv = CSLAddString(argv, "-r"); // resampling
899  argv = CSLAddString(argv, interp.c_str()); // {nearest(default),bilinear,cubic,cubicspline,lanczos,average,mode}
900 
901  if (!crs.empty()) {
902  argv = CSLAddString(argv, "-a_srs"); // dst SRS (WKT or "EPSG:n")
903  argv = CSLAddString(argv, crs.c_str());
904  }
905 
906  if (BESISDEBUG(DEBUG_KEY)) {
907  char **local = argv;
908  while (*local) {
909  BESDEBUG(DEBUG_KEY, "argv: " << *local++ << endl);
910  }
911 
912  }
913 #if 0
914  char **local = argv;
915  while (*local) {
916  cerr << "argv: " << *local++ << endl;
917  }
918 #endif
919 
920  GDALTranslateOptions *options = GDALTranslateOptionsNew(argv, NULL /*binary options*/);
921 
922  int usage_error = CE_None; // result
923  GDALDatasetH dst_handle = GDALTranslate("warped_dst", src.get(), options, &usage_error);
924  if (!dst_handle || usage_error != CE_None) {
925  GDALClose(dst_handle);
926  GDALTranslateOptionsFree(options);
927  string msg = string("Error calling GDAL translate: ") + CPLGetLastErrorMsg();
928  BESDEBUG(DEBUG_KEY, "ERROR scale_dataset(): " << msg << endl);
929  throw BESError(msg, BES_INTERNAL_ERROR, __FILE__, __LINE__);
930  }
931 
932  unique_ptr<GDALDataset> dst(static_cast<GDALDataset*>(dst_handle));
933 
934  GDALTranslateOptionsFree(options);
935 
936  return dst;
937 }
938 
939 
940 unique_ptr<GDALDataset> scale_dataset_3D(unique_ptr<GDALDataset>& src, const SizeBox &size, const string &crs /*""*/,
941  const string &interp /*nearest*/)
942 {
943  char **argv = NULL;
944  argv = CSLAddString(argv, "-of"); // output format
945  argv = CSLAddString(argv, "MEM");
946 
947  argv = CSLAddString(argv, "-outsize"); // output size
948  ostringstream oss;
949  oss << size.x_size;
950  argv = CSLAddString(argv, oss.str().c_str()); // size x
951  oss.str("");
952  oss << size.y_size;
953  argv = CSLAddString(argv, oss.str().c_str()); // size y
954 
955  // all bands in src
956  int n_bands = src.get()->GetRasterCount();
957  for(int i=0; i < n_bands; i++){
958  oss.str("");
959  oss << i+1;
960  argv = CSLAddString(argv, "-b");
961  argv = CSLAddString(argv, oss.str().c_str());
962  }
963 
964  argv = CSLAddString(argv, "-r"); // resampling
965  argv = CSLAddString(argv, interp.c_str()); // {nearest(default),bilinear,cubic,cubicspline,lanczos,average,mode}
966 
967  if (!crs.empty()) {
968  argv = CSLAddString(argv, "-a_srs"); // dst SRS (WKT or "EPSG:n")
969  argv = CSLAddString(argv, crs.c_str());
970  }
971 
972  if (BESISDEBUG(DEBUG_KEY)) {
973  char **local = argv;
974  while (*local) {
975  BESDEBUG(DEBUG_KEY, "argv: " << *local++ << endl);
976  }
977 
978  }
979 #if 0
980  char **local = argv;
981  while (*local) {
982  cerr << "argv: " << *local++ << endl;
983  }
984 #endif
985 
986  GDALTranslateOptions *options = GDALTranslateOptionsNew(argv, NULL /*binary options*/);
987  int usage_error = CE_None; // result
988  GDALDatasetH dst_handle = GDALTranslate("warped_dst", src.get(), options, &usage_error);
989  if (!dst_handle || usage_error != CE_None) {
990  GDALClose(dst_handle);
991  GDALTranslateOptionsFree(options);
992  string msg = string("Error calling GDAL translate: ") + CPLGetLastErrorMsg();
993  BESDEBUG(DEBUG_KEY, "ERROR scale_dataset(): " << msg << endl);
994  throw BESError(msg, BES_INTERNAL_ERROR, __FILE__, __LINE__);
995  }
996 
997  unique_ptr<GDALDataset> dst(static_cast<GDALDataset*>(dst_handle));
998 
999  GDALTranslateOptionsFree(options);
1000 
1001  return dst;
1002 }
1003 
1004 
1017 Grid *scale_dap_array(const Array *data, const Array *x, const Array *y, const SizeBox &size,
1018  const string &crs, const string &interp)
1019 {
1020  // Build GDALDataset for Grid g with lon and lat maps as given
1021  Array *d = const_cast<Array*>(data);
1022 
1023  unique_ptr<GDALDataset> src = build_src_dataset(d, const_cast<Array*>(x), const_cast<Array*>(y));
1024 
1025  // scale to the new size, using optional CRS and interpolation params
1026  unique_ptr<GDALDataset> dst = scale_dataset(src, size, crs, interp);
1027 
1028  // Build a result Grid: extract the data, build the maps and assemble
1029  unique_ptr<Array> built_data(build_array_from_gdal_dataset(dst.get(), d));
1030 
1031  unique_ptr<Array> built_lat(new Array(y->name(), new Float32(y->name())));
1032  unique_ptr<Array> built_lon(new Array(x->name(), new Float32(x->name())));
1033 
1034  build_maps_from_gdal_dataset(dst.get(), built_lon.get(), built_lat.get());
1035 
1036  unique_ptr<Grid> result(new Grid(d->name()));
1037  result->set_array(built_data.release());
1038  result->add_map(built_lat.release(), false);
1039  result->add_map(built_lon.release(), false);
1040 
1041  return result.release();
1042 }
1043 
1054 Grid *scale_dap_grid(const Grid *g, const SizeBox &size, const string &crs, const string &interp)
1055 {
1056  string func(__func__);
1057  func += "() - ";
1058 
1059  if(!g){
1060  throw BESError(func+"The Grid object is null.",BES_INTERNAL_ERROR,__FILE__,__LINE__);
1061  }
1062  // Build GDALDataset for Grid g with lon and lat maps as given
1063  Array *data = dynamic_cast<Array*>(const_cast<Grid*>(g)->array_var());
1064  if(!data){
1065  throw BESError(func+"Unable to obtain data array from Grid '"+g->name()+"'",BES_INTERNAL_ERROR,__FILE__,__LINE__);
1066  }
1067  // return iteration
1068  Grid::Map_riter ritr = const_cast<Grid*>(g)->map_rbegin();
1069  Array *x = dynamic_cast<Array*>(*ritr);
1070  Array *y = dynamic_cast<Array*>(*++ritr);
1071 
1072  if(!x || !y){
1073  throw BESError(func+"Unable to obtain 2 Map arrays from Grid '"+g->name()+"'",BES_INTERNAL_ERROR,__FILE__,__LINE__);
1074  }
1075 
1076  return scale_dap_array(data, x, y, size, crs, interp);
1077 }
1078 
1079 #define ADD_BAND 0
1080 
1099 unique_ptr<GDALDataset> build_src_dataset_3D(Array *data, Array *t, Array *x, Array *y, const string &srs)
1100 {
1101  Array *d = dynamic_cast<Array*>(data);
1102 
1103  GDALDriver *driver = GetGDALDriverManager()->GetDriverByName("MEM");
1104  if(!driver){
1105  string msg = string("Could not get the Memory driver for GDAL: ") + CPLGetLastErrorMsg();
1106  BESDEBUG(DEBUG_KEY, "ERROR build_src_dataset(): " << msg << endl);
1107  throw BESError(msg,BES_INTERNAL_ERROR,__FILE__,__LINE__);
1108  }
1109 
1110  SizeBox array_size = get_size_box(x, y);
1111  int nBands = t->length();
1112  BESDEBUG(DEBUG_KEY, "nBands = " << nBands << endl);
1113  int nBytes = data->prototype()->width();
1114  const int data_size = x->length() * y->length();
1115  unsigned int dsize = data_size * nBytes;
1116 
1117  unique_ptr<GDALDataset> ds(driver->Create("result", array_size.x_size, array_size.y_size, nBands, get_array_type(d),
1118  NULL /* driver_options */));
1119  data->read();
1120  // start band loop
1121  for(int i=1; i<=nBands; i++){
1122 
1123  GDALRasterBand *band = ds->GetRasterBand(i);
1124  if (!band) {
1125  string msg = "Could not get the GDAL RasterBand for Array '" + data->name() + "': " + CPLGetLastErrorMsg();
1126  BESDEBUG(DEBUG_KEY,"ERROR build_src_dataset(): " << msg << endl);
1127  throw BESError(msg,BES_INTERNAL_ERROR,__FILE__,__LINE__);
1128  }
1129 
1130  double no_data = get_missing_data_value(data);
1131  band->SetNoDataValue(no_data);
1132 
1133  #if !ADD_BAND
1134  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);
1135  if (error != CPLE_None)
1136  throw Error("Could not write data for band: " + long_to_string(i) + ": " + string(CPLGetLastErrorMsg()));
1137  #endif
1138 
1139 
1140  } // end band loop
1141  vector<double> geo_transform = get_geotransform_data(x, y);
1142  ds->SetGeoTransform(&geo_transform[0]);
1143 
1144  OGRSpatialReference native_srs;
1145  if (CE_None != native_srs.SetWellKnownGeogCS(srs.c_str())){
1146  string msg = "Could not set '" + srs + "' as the dataset native CRS.";
1147  BESDEBUG(DEBUG_KEY,"ERROR build_src_dataset(): " << msg << endl);
1148  throw BESError(msg,BES_SYNTAX_USER_ERROR,__FILE__,__LINE__);
1149  }
1150 
1151  // Connect the SRS/CRS to the GDAL Dataset
1152  char *pszSRS_WKT = NULL;
1153  native_srs.exportToWkt( &pszSRS_WKT );
1154  ds->SetProjection( pszSRS_WKT );
1155  CPLFree( pszSRS_WKT );
1156 
1157  return ds;
1158 }
1159 
1173 Grid *scale_dap_array_3D(const Array *data, const Array *time, const Array *lon, const Array *lat, const SizeBox &size,
1174  const string &crs, const string &interp)
1175 {
1176  // Build GDALDataset for Grid g with lon and lat maps as given
1177  Array *d = const_cast<Array*>(data);
1178  Array *t = const_cast<Array*>(time);
1179  Array *x = const_cast<Array*>(lon);
1180  Array *y = const_cast<Array*>(lat);
1181 
1182  // get GDALDataset with bands
1183  unique_ptr<GDALDataset> src = build_src_dataset_3D(d, const_cast<Array*>(t), const_cast<Array*>(x), const_cast<Array*>(y));
1184 
1185  // scale all bands to the new size, using optional CRS and interpolation params
1186  unique_ptr<GDALDataset> dst = scale_dataset_3D(src, size, crs, interp);
1187 
1188  // Build a result Grid: extract the data, build the maps and assemble
1189  unique_ptr<Array> built_data(build_array_from_gdal_dataset_3D(dst.get(), d));
1190  unique_ptr<Array> built_time(new Array(t->name(), new Float32(t->name())));
1191  unique_ptr<Array> built_lat(new Array(y->name(), new Float32(y->name())));
1192  unique_ptr<Array> built_lon(new Array(x->name(), new Float32(x->name())));
1193 
1194  // Build maps for grid
1195  build_maps_from_gdal_dataset_3D(dst.get(), t, built_time.get(), built_lon.get(), built_lat.get());
1196 
1197  // get result Grid
1198  unique_ptr<Grid> result(new Grid(d->name()));
1199  result->set_array(built_data.release());
1200  result->add_map(built_time.release(), false);
1201  result->add_map(built_lat.release(), false);
1202  result->add_map(built_lon.release(), false);
1203 
1204  return result.release();
1205 }
1206 
1207 }
1208 
libdap
Definition: BESDapFunctionResponseCache.h:35
Error
BESError
Abstract exception class for the BES with basic string message.
Definition: BESError.h:58