bes  Updated for version 3.20.6
gdal_utils.cc
1 // This file is part of the GDAL OPeNDAP Adapter
2 
3 // Copyright (c) 2004 OPeNDAP, Inc.
4 // Author: Frank Warmerdam <warmerdam@pobox.com>
5 //
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License, or (at your option) any later version.
10 //
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
15 //
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
19 //
20 // You can contact OPeNDAP, Inc. at PO Box 112, Saunderstown, RI. 02874-0112.
21 
22 #include "config.h"
23 
24 #include <iostream>
25 #include <sstream>
26 #include <string>
27 
28 #include <gdal.h>
29 #include <cpl_string.h>
30 
31 //#define DODS_DEBUG 1
32 
33 #include <Byte.h>
34 #include <UInt16.h>
35 #include <Int16.h>
36 #include <UInt32.h>
37 #include <Int32.h>
38 #include <Float32.h>
39 #include <Float64.h>
40 
41 #include <DDS.h>
42 #include <DAS.h>
43 #include <BaseTypeFactory.h>
44 
45 #include <DMR.h>
46 #include <D4Group.h>
47 #include <D4Dimensions.h>
48 #include <D4Maps.h>
49 #include <D4Attributes.h>
50 #include <D4BaseTypeFactory.h>
51 #include <debug.h>
52 
53 #include "GDALTypes.h"
54 
55 using namespace libdap;
56 
57 /************************************************************************/
58 /* attach_str_attr_item() */
59 /* */
60 /* Add a string attribute item to target container with */
61 /* appropriate quoting and escaping. */
62 /************************************************************************/
63 
64 static void attach_str_attr_item(AttrTable *parent_table, const char *pszKey, const char *pszValue)
65 {
66  //string oQuotedValue;
67  char *pszEscapedText = CPLEscapeString(pszValue, -1,
68  CPLES_BackslashQuotable);
69 
70  parent_table->append_attr(pszKey, "String", pszEscapedText /*oQuotedValue*/);
71 
72  CPLFree(pszEscapedText);
73 }
74 
75 /************************************************************************/
76 /* translate_metadata() */
77 /* */
78 /* Turn a list of metadata name/value pairs into DAS into and */
79 /* attach it to the passed container. */
80 /************************************************************************/
81 
82 static void translate_metadata(char **md, AttrTable *parent_table)
83 {
84  AttrTable *md_table;
85  int i;
86 
87  md_table = parent_table->append_container(string("Metadata"));
88 
89  for (i = 0; md != NULL && md[i] != NULL; i++) {
90  const char *pszValue;
91  char *pszKey = NULL;
92 
93  pszValue = CPLParseNameValue(md[i], &pszKey);
94 
95  attach_str_attr_item(md_table, pszKey, pszValue);
96 
97  CPLFree(pszKey);
98  }
99 }
100 
107 static void build_global_attributes(const GDALDatasetH& hDS, AttrTable* attr_table)
108 {
109  /* -------------------------------------------------------------------- */
110  /* Geotransform */
111  /* -------------------------------------------------------------------- */
112  double adfGeoTransform[6];
113  if (GDALGetGeoTransform(hDS, adfGeoTransform) == CE_None
114  && (adfGeoTransform[0] != 0.0 || adfGeoTransform[1] != 1.0 || adfGeoTransform[2] != 0.0
115  || adfGeoTransform[3] != 0.0 || adfGeoTransform[4] != 0.0 || fabs(adfGeoTransform[5]) != 1.0)) {
116 
117  char szGeoTransform[200];
118  double dfMaxX, dfMinX, dfMaxY, dfMinY;
119  int nXSize = GDALGetRasterXSize(hDS);
120  int nYSize = GDALGetRasterYSize(hDS);
121 
122  dfMaxX =
123  MAX(MAX(adfGeoTransform[0], adfGeoTransform[0] + adfGeoTransform[1] * nXSize),
124  MAX(adfGeoTransform[0] + adfGeoTransform[2] * nYSize, adfGeoTransform[0] + adfGeoTransform[2] * nYSize + adfGeoTransform[1] * nXSize));
125  dfMinX =
126  MIN(MIN(adfGeoTransform[0], adfGeoTransform[0] + adfGeoTransform[1] * nXSize),
127  MIN(adfGeoTransform[0] + adfGeoTransform[2] * nYSize, adfGeoTransform[0] + adfGeoTransform[2] * nYSize + adfGeoTransform[1] * nXSize));
128  dfMaxY =
129  MAX(MAX(adfGeoTransform[3], adfGeoTransform[3] + adfGeoTransform[4] * nXSize),
130  MAX(adfGeoTransform[3] + adfGeoTransform[5] * nYSize, adfGeoTransform[3] + adfGeoTransform[5] * nYSize + adfGeoTransform[4] * nXSize));
131  dfMinY =
132  MIN(MIN(adfGeoTransform[3], adfGeoTransform[3] + adfGeoTransform[4] * nXSize),
133  MIN(adfGeoTransform[3] + adfGeoTransform[5] * nYSize, adfGeoTransform[3] + adfGeoTransform[5] * nYSize + adfGeoTransform[4] * nXSize));
134 
135  attr_table->append_attr("Northernmost_Northing", "Float64", CPLSPrintf("%.16g", dfMaxY));
136  attr_table->append_attr("Southernmost_Northing", "Float64", CPLSPrintf("%.16g", dfMinY));
137  attr_table->append_attr("Easternmost_Easting", "Float64", CPLSPrintf("%.16g", dfMaxX));
138 #if 0
139  // Gareth Williams pointed out this typo. jhrg 9/26/19
140  attr_table->append_attr("Westernmost_Northing", "Float64", CPLSPrintf("%.16g", dfMinX));
141 #endif
142  attr_table->append_attr("Westernmost_Easting", "Float64", CPLSPrintf("%.16g", dfMinX));
143 
144  snprintf(szGeoTransform, 200, "%.16g %.16g %.16g %.16g %.16g %.16g", adfGeoTransform[0], adfGeoTransform[1],
145  adfGeoTransform[2], adfGeoTransform[3], adfGeoTransform[4], adfGeoTransform[5]);
146 
147  attach_str_attr_item(attr_table, "GeoTransform", szGeoTransform);
148  }
149 
150  /* -------------------------------------------------------------------- */
151  /* Metadata. */
152  /* -------------------------------------------------------------------- */
153  char** md;
154  md = GDALGetMetadata(hDS, NULL);
155  if (md != NULL) translate_metadata(md, attr_table);
156 
157  /* -------------------------------------------------------------------- */
158  /* SRS */
159  /* -------------------------------------------------------------------- */
160  const char* pszWKT = GDALGetProjectionRef(hDS);
161  if (pszWKT != NULL && strlen(pszWKT) > 0) attach_str_attr_item(attr_table, "spatial_ref", pszWKT);
162 }
163 
171 static void build_variable_attributes(const GDALDatasetH &hDS, AttrTable *band_attr, const int iBand)
172 {
173  /* -------------------------------------------------------------------- */
174  /* Offset. */
175  /* -------------------------------------------------------------------- */
176  int bSuccess;
177  double dfValue;
178  char szValue[128];
179  GDALRasterBandH hBand = GDALGetRasterBand(hDS, iBand + 1);
180 
181  dfValue = GDALGetRasterOffset(hBand, &bSuccess);
182  if (bSuccess) {
183  snprintf(szValue, 128, "%.16g", dfValue);
184  band_attr->append_attr("add_offset", "Float64", szValue);
185  }
186 
187  /* -------------------------------------------------------------------- */
188  /* Scale */
189  /* -------------------------------------------------------------------- */
190  dfValue = GDALGetRasterScale(hBand, &bSuccess);
191  if (bSuccess) {
192  snprintf(szValue, 128, "%.16g", dfValue);
193  band_attr->append_attr("scale_factor", "Float64", szValue);
194  }
195 
196  /* -------------------------------------------------------------------- */
197  /* nodata/missing_value */
198  /* -------------------------------------------------------------------- */
199  dfValue = GDALGetRasterNoDataValue(hBand, &bSuccess);
200  if (bSuccess) {
201  snprintf(szValue, 128, "%.16g", dfValue);
202  band_attr->append_attr("missing_value", "Float64", szValue);
203  }
204 
205  /* -------------------------------------------------------------------- */
206  /* Description. */
207  /* -------------------------------------------------------------------- */
208  if (GDALGetDescription(hBand) != NULL && strlen(GDALGetDescription(hBand)) > 0) {
209  attach_str_attr_item(band_attr, "Description", GDALGetDescription(hBand));
210  }
211 
212  /* -------------------------------------------------------------------- */
213  /* PhotometricInterpretation. */
214  /* -------------------------------------------------------------------- */
215  if (GDALGetRasterColorInterpretation(hBand) != GCI_Undefined) {
216  attach_str_attr_item(band_attr, "PhotometricInterpretation",
217  GDALGetColorInterpretationName(GDALGetRasterColorInterpretation(hBand)));
218  }
219 
220  /* -------------------------------------------------------------------- */
221  /* Band Metadata. */
222  /* -------------------------------------------------------------------- */
223  char **md = GDALGetMetadata(hBand, NULL);
224  if (md != NULL) translate_metadata(md, band_attr);
225 
226  /* -------------------------------------------------------------------- */
227  /* Colormap. */
228  /* -------------------------------------------------------------------- */
229  GDALColorTableH hCT;
230 
231  hCT = GDALGetRasterColorTable(hBand);
232  if (hCT != NULL) {
233  AttrTable *ct_attr;
234  int iColor, nColorCount = GDALGetColorEntryCount(hCT);
235 
236  ct_attr = band_attr->append_container(string("Colormap"));
237 
238  for (iColor = 0; iColor < nColorCount; iColor++) {
239  GDALColorEntry sRGB;
240  AttrTable *color_attr;
241 
242  color_attr = ct_attr->append_container(string(CPLSPrintf("color_%d", iColor)));
243 
244  GDALGetColorEntryAsRGB(hCT, iColor, &sRGB);
245 
246  color_attr->append_attr("red", "byte", CPLSPrintf("%d", sRGB.c1));
247  color_attr->append_attr("green", "byte", CPLSPrintf("%d", sRGB.c2));
248  color_attr->append_attr("blue", "byte", CPLSPrintf("%d", sRGB.c3));
249  color_attr->append_attr("alpha", "byte", CPLSPrintf("%d", sRGB.c4));
250  }
251  }
252 }
253 
267 void gdal_read_dataset_attributes(DAS &das, const GDALDatasetH &hDS)
268 {
269  AttrTable *attr_table = das.add_table(string("GLOBAL"), new AttrTable);
270 
271  build_global_attributes(hDS, attr_table);
272 
273  /* ==================================================================== */
274  /* Generation info for bands. */
275  /* ==================================================================== */
276  for (int iBand = 0; iBand < GDALGetRasterCount(hDS); iBand++) {
277  char szName[128];
278  AttrTable *band_attr;
279 
280  /* -------------------------------------------------------------------- */
281  /* Create container named after the band. */
282  /* -------------------------------------------------------------------- */
283  snprintf(szName, 128, "band_%d", iBand + 1);
284  band_attr = das.add_table(string(szName), new AttrTable);
285 
286  build_variable_attributes(hDS, band_attr, iBand);
287  }
288 }
289 
300 void gdal_read_dataset_variables(DDS *dds, const GDALDatasetH &hDS, const string &filename,bool include_attrs)
301 {
302  // Load in to global attributes
303  if(true == include_attrs) {
304  AttrTable *global_attr = dds->get_attr_table().append_container("GLOBAL");
305  build_global_attributes(hDS, global_attr);
306  }
307 
308  /* -------------------------------------------------------------------- */
309  /* Create the basic matrix for each band. */
310  /* -------------------------------------------------------------------- */
311  BaseTypeFactory factory; // Use this to build new scalar variables
312  GDALDataType eBufType;
313 
314  for (int iBand = 0; iBand < GDALGetRasterCount(hDS); iBand++) {
315  DBG(cerr << "In dgal_dds.cc iBand" << endl);
316 
317  GDALRasterBandH hBand = GDALGetRasterBand(hDS, iBand + 1);
318 
319  ostringstream oss;
320  oss << "band_" << iBand + 1;
321 
322  eBufType = GDALGetRasterDataType(hBand);
323 
324  BaseType *bt;
325  switch (GDALGetRasterDataType(hBand)) {
326  case GDT_Byte:
327  bt = factory.NewByte(oss.str());
328  break;
329 
330  case GDT_UInt16:
331  bt = factory.NewUInt16(oss.str());
332  break;
333 
334  case GDT_Int16:
335  bt = factory.NewInt16(oss.str());
336  break;
337 
338  case GDT_UInt32:
339  bt = factory.NewUInt32(oss.str());
340  break;
341 
342  case GDT_Int32:
343  bt = factory.NewInt32(oss.str());
344  break;
345 
346  case GDT_Float32:
347  bt = factory.NewFloat32(oss.str());
348  break;
349 
350  case GDT_Float64:
351  bt = factory.NewFloat64(oss.str());
352  break;
353 
354  case GDT_CFloat32:
355  case GDT_CFloat64:
356  case GDT_CInt16:
357  case GDT_CInt32:
358  default:
359  // TODO eventually we need to preserve complex info
360  bt = factory.NewFloat64(oss.str());
361  eBufType = GDT_Float64;
362  break;
363  }
364 
365  /* -------------------------------------------------------------------- */
366  /* Create a grid to hold the raster. */
367  /* -------------------------------------------------------------------- */
368  Grid *grid;
369  grid = new GDALGrid(filename, oss.str());
370 
371  /* -------------------------------------------------------------------- */
372  /* Make into an Array for the raster data with appropriate */
373  /* dimensions. */
374  /* -------------------------------------------------------------------- */
375  Array *ar;
376  // A 'feature' of Array is that it copies the variable passed to
377  // its ctor. To get around that, pass null and use add_var_nocopy().
378  // Modified for the DAP4 response; switched to this new ctor.
379  ar = new GDALArray(oss.str(), 0, filename, eBufType, iBand + 1);
380  ar->add_var_nocopy(bt);
381  ar->append_dim(GDALGetRasterYSize(hDS), "northing");
382  ar->append_dim(GDALGetRasterXSize(hDS), "easting");
383 
384  grid->add_var_nocopy(ar, libdap::array);
385 
386  /* -------------------------------------------------------------------- */
387  /* Add the dimension map arrays. */
388  /* -------------------------------------------------------------------- */
389  //bt = new GDALFloat64( "northing" );
390  bt = factory.NewFloat64("northing");
391  ar = new GDALArray("northing", 0, filename, GDT_Float64/* eBufType */, iBand + 1);
392  ar->add_var_nocopy(bt);
393  ar->append_dim(GDALGetRasterYSize(hDS), "northing");
394 
395  grid->add_var_nocopy(ar, maps);
396 
397  //bt = new GDALFloat64( "easting" );
398  bt = factory.NewFloat64("easting");
399  ar = new GDALArray("easting", 0, filename, GDT_Float64/* eBufType */, iBand + 1);
400  ar->add_var_nocopy(bt);
401  ar->append_dim(GDALGetRasterXSize(hDS), "easting");
402 
403  grid->add_var_nocopy(ar, maps);
404 
405  DBG(cerr << "Type of grid: " << typeid(grid).name() << endl);
406 
407  // Add attributes to the Grid
408  if(true == include_attrs) {
409  AttrTable &band_attr = grid->get_attr_table();
410  build_variable_attributes(hDS, &band_attr, iBand);
411  }
412 
413  dds->add_var_nocopy(grid);
414  }
415 }
416 
430 void gdal_read_dataset_variables(DMR *dmr, const GDALDatasetH &hDS, const string &filename)
431 {
432  // Load in global attributes. Hack, use DAP2 attributes and transform.
433  // This is easier than writing new D4Attributes code and not a heuristic
434  // routine like the variable transforms or attribute to DDS merge. New
435  // code for the D4Attributes would duplicate the AttrTable code.
436  //
437  // An extra AttrTable is needed because of oddities of its API but that
438  // comes in handy since D4Attributes::transform_to_dap4() throws away
439  // the top most AttrTable
440 
441  AttrTable *attr = new AttrTable();
442  AttrTable *global_attr = attr->append_container("GLOBAL");
443 
444  build_global_attributes(hDS, global_attr);
445 
446  dmr->root()->attributes()->transform_to_dap4(*attr);
447  delete attr; attr = 0;
448 
449  D4BaseTypeFactory factory; // Use this to build new variables
450 
451  // Make the northing and easting shared dims for this gdal dataset.
452  // For bands that have a different size than the overall Raster{X,Y}Size,
453  // assume they are lower resolution bands. For these I'm going to simply
454  // not include the shared dimensions for them. If this is a problem,
455  // we can expand the set of dimensions later.
456  //
457  // See the GDAL data model doc (http://www.gdal.org/gdal_datamodel.html)
458  // for info on why this seems correct. jhrg 6/2/17
459 
460  // Find the first band that is the same size as the whole raster dataset (i.e.,
461  // the first band that is not at a reduced resolution). In the larger loop
462  // below we only bind sdims to bands that match the overall raster size and
463  // leave bands that are at a reduce resolution w/o shared dims.
464  int sdim_band_num = 1;
465  for (; sdim_band_num <= GDALGetRasterCount(hDS); ++sdim_band_num) {
466  if (GDALGetRasterBandYSize(GDALGetRasterBand(hDS, sdim_band_num)) == GDALGetRasterYSize(hDS)
467  && GDALGetRasterBandXSize(GDALGetRasterBand(hDS, sdim_band_num)) == GDALGetRasterXSize(hDS)) {
468  break;
469  }
470  }
471 
472  // Make the two shared dims that will have the geo-referencing info
473  const string northing = "northing";
474  // Add the shared dim
475  D4Dimension *northing_dim = new D4Dimension(northing, GDALGetRasterYSize(hDS));
476  dmr->root()->dims()->add_dim_nocopy(northing_dim);
477  // use the shared dim for the map
478  Array *northing_map = new GDALArray(northing, 0, filename, GDT_Float64, sdim_band_num);
479  northing_map->add_var_nocopy(factory.NewFloat64(northing));
480  northing_map->append_dim(northing_dim);
481  // add the map
482  dmr->root()->add_var_nocopy(northing_map); // Add the map array to the DMR/D4Group
483 
484  const string easting = "easting";
485  D4Dimension *easting_dim = new D4Dimension(easting, GDALGetRasterXSize(hDS));
486  dmr->root()->dims()->add_dim_nocopy(easting_dim);
487  Array *easting_map = new GDALArray(easting, 0, filename, GDT_Float64, sdim_band_num);
488  easting_map->add_var_nocopy(factory.NewFloat64(easting));
489  easting_map->append_dim(easting_dim);
490  dmr->root()->add_var_nocopy(easting_map); // Add the map array to the DMR/D4Group
491 
492  /* -------------------------------------------------------------------- */
493  /* Create the basic matrix for each band. */
494  /* -------------------------------------------------------------------- */
495 
496  GDALDataType eBufType;
497 
498  for (int iBand = 0; iBand < GDALGetRasterCount(hDS); iBand++) {
499  DBG(cerr << "In dgal_dds.cc iBand" << endl);
500 
501  GDALRasterBandH hBand = GDALGetRasterBand(hDS, iBand + 1);
502 
503  ostringstream oss;
504  oss << "band_" << iBand + 1;
505 
506  eBufType = GDALGetRasterDataType(hBand);
507 
508  BaseType *bt;
509  switch (GDALGetRasterDataType(hBand)) {
510  case GDT_Byte:
511  bt = factory.NewByte(oss.str());
512  break;
513 
514  case GDT_UInt16:
515  bt = factory.NewUInt16(oss.str());
516  break;
517 
518  case GDT_Int16:
519  bt = factory.NewInt16(oss.str());
520  break;
521 
522  case GDT_UInt32:
523  bt = factory.NewUInt32(oss.str());
524  break;
525 
526  case GDT_Int32:
527  bt = factory.NewInt32(oss.str());
528  break;
529 
530  case GDT_Float32:
531  bt = factory.NewFloat32(oss.str());
532  break;
533 
534  case GDT_Float64:
535  bt = factory.NewFloat64(oss.str());
536  break;
537 
538  case GDT_CFloat32:
539  case GDT_CFloat64:
540  case GDT_CInt16:
541  case GDT_CInt32:
542  default:
543  // TODO eventually we need to preserve complex info
544  // Replace this with an attribute about a missing/elided variable?
545  bt = factory.NewFloat64(oss.str());
546  eBufType = GDT_Float64;
547  break;
548  }
549 
550  // Make the array for this band and then associate the dimensions/maps
551  // once they are made;
552  Array *ar = new GDALArray(oss.str(), 0, filename, eBufType, iBand + 1);
553  ar->add_var_nocopy(bt); // bt is from the above switch
554 
555  /* -------------------------------------------------------------------- */
556  /* Add the dimension map arrays. */
557  /* -------------------------------------------------------------------- */
558 
559  if (GDALGetRasterBandYSize(hBand) == GDALGetRasterYSize(hDS)
560  && GDALGetRasterBandXSize(hBand) == GDALGetRasterXSize(hDS)) {
561  // Use the shared dimensions
562  ar->append_dim(northing_dim);
563  ar->append_dim(easting_dim);
564 
565  // Bind the map to the array; FQNs for the maps (shared dims)
566  ar->maps()->add_map(new D4Map(string("/") + northing, northing_map, ar));
567  ar->maps()->add_map(new D4Map(string("/") + easting, easting_map, ar));
568  }
569  else {
570  // Don't use the shared dims
571  ar->append_dim(GDALGetRasterBandYSize(hBand), northing);
572  ar->append_dim(GDALGetRasterBandXSize(hBand), easting);
573  }
574 
575  // Add attributes, using the same hack as for the global attributes;
576  // build the DAP2 AttrTable object and then transform to DAP4.
577  AttrTable &band_attr = ar->get_attr_table();
578  build_variable_attributes(hDS, &band_attr, iBand);
579  ar->attributes()->transform_to_dap4(band_attr);
580 
581  dmr->root()->add_var_nocopy(ar); // Add the array to the DMR
582  }
583 }
584 
602 void read_data_array(GDALArray *array, const GDALRasterBandH &hBand)
603 {
604  /* -------------------------------------------------------------------- */
605  /* Collect the x and y sampling values from the constraint. */
606  /* -------------------------------------------------------------------- */
607  Array::Dim_iter p = array->dim_begin();
608  int start = array->dimension_start(p, true);
609  int stride = array->dimension_stride(p, true);
610  int stop = array->dimension_stop(p, true);
611 
612  // Test for the case where a dimension has not been subset. jhrg 2/18/16
613  if (array->dimension_size(p, true) == 0) { //default rows
614  start = 0;
615  stride = 1;
616  stop = GDALGetRasterBandYSize(hBand) - 1;
617  }
618 
619  p++;
620  int start_2 = array->dimension_start(p, true);
621  int stride_2 = array->dimension_stride(p, true);
622  int stop_2 = array->dimension_stop(p, true);
623 
624  if (array->dimension_size(p, true) == 0) { //default columns
625  start_2 = 0;
626  stride_2 = 1;
627  stop_2 = GDALGetRasterBandXSize(hBand) - 1;
628  }
629 
630  /* -------------------------------------------------------------------- */
631  /* Build a window and buf size from this. */
632  /* -------------------------------------------------------------------- */
633  int nWinXOff, nWinYOff, nWinXSize, nWinYSize, nBufXSize, nBufYSize;
634 
635  nWinXOff = start_2;
636  nWinYOff = start;
637  nWinXSize = stop_2 + 1 - start_2;
638  nWinYSize = stop + 1 - start;
639 
640  nBufXSize = (stop_2 - start_2) / stride_2 + 1;
641  nBufYSize = (stop - start) / stride + 1;
642 
643  /* -------------------------------------------------------------------- */
644  /* Allocate buffer. */
645  /* -------------------------------------------------------------------- */
646  int nPixelSize = GDALGetDataTypeSize(array->get_gdal_buf_type()) / 8;
647  vector<char> pData(nBufXSize * nBufYSize * nPixelSize);
648 
649  /* -------------------------------------------------------------------- */
650  /* Read request into buffer. */
651  /* -------------------------------------------------------------------- */
652  CPLErr eErr = GDALRasterIO(hBand, GF_Read, nWinXOff, nWinYOff, nWinXSize, nWinYSize, &pData[0], nBufXSize,
653  nBufYSize, array->get_gdal_buf_type(), 0, 0);
654  if (eErr != CE_None) throw Error("Error reading: " + array->name());
655 
656  array->val2buf(&pData[0]);
657 }
658 
677 void read_map_array(Array *map, const GDALRasterBandH &hBand, const GDALDatasetH &hDS)
678 {
679  Array::Dim_iter p = map->dim_begin();
680  int start = map->dimension_start(p, true);
681  int stride = map->dimension_stride(p, true);
682  int stop = map->dimension_stop(p, true);
683 
684  if (start + stop + stride == 0) { //default rows
685  start = 0;
686  stride = 1;
687  if (map->name() == "northing")
688  stop = GDALGetRasterBandYSize(hBand) - 1;
689  else if (map->name() == "easting")
690  stop = GDALGetRasterBandXSize(hBand) - 1;
691  else
692  throw Error("Expected a map named 'northing' or 'easting' but got: " + map->name());
693  }
694 
695  int nBufSize = (stop - start) / stride + 1;
696 
697  /* -------------------------------------------------------------------- */
698  /* Read or default the geotransform used to generate the */
699  /* georeferencing maps. */
700  /* -------------------------------------------------------------------- */
701 
702  double adfGeoTransform[6];
703 
704  if (GDALGetGeoTransform(hDS, adfGeoTransform) != CE_None) {
705  adfGeoTransform[0] = 0.0;
706  adfGeoTransform[1] = 1.0;
707  adfGeoTransform[2] = 0.0;
708  adfGeoTransform[3] = 0.0;
709  adfGeoTransform[4] = 0.0;
710  adfGeoTransform[5] = 1.0;
711  }
712 
713  /* -------------------------------------------------------------------- */
714  /* Set the map array. */
715  /* -------------------------------------------------------------------- */
716  vector<double> padfMap(nBufSize);
717 
718  if (map->name() == "northing") {
719  for (int i = 0, iLine = start; iLine <= stop; iLine += stride) {
720  padfMap[i++] = adfGeoTransform[3] + adfGeoTransform[5] * iLine;
721  }
722  }
723  else if (map->name() == "easting") {
724  for (int i = 0, iPixel = start; iPixel <= stop; iPixel += stride) {
725  padfMap[i++] = adfGeoTransform[0] + iPixel * adfGeoTransform[1];
726  }
727  }
728  else
729  throw Error("Expected a map named 'northing' or 'easting' but got: " + map->name());
730 
731  map->val2buf((void *) &padfMap[0]);
732 }
733 
734 // $Log: gdal_das.cc,v $
735 // Revision 1.1 2004/10/19 20:38:28 warmerda
736 // New
737 //
738 // Revision 1.2 2004/10/15 18:06:45 warmerda
739 // Strip the extension off the filename.
740 //
741 // Revision 1.1 2004/10/04 14:29:29 warmerda
742 // New
743 //
744 
libdap
Definition: BESDapFunctionResponseCache.h:35
GDALGrid
Definition: GDALTypes.h:59
GDALArray
Definition: GDALTypes.h:34
Error