bes  Updated for version 3.20.5
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  attr_table->append_attr("Westernmost_Northing", "Float64", CPLSPrintf("%.16g", dfMinX));
139 
140  snprintf(szGeoTransform, 200, "%.16g %.16g %.16g %.16g %.16g %.16g", adfGeoTransform[0], adfGeoTransform[1],
141  adfGeoTransform[2], adfGeoTransform[3], adfGeoTransform[4], adfGeoTransform[5]);
142 
143  attach_str_attr_item(attr_table, "GeoTransform", szGeoTransform);
144  }
145 
146  /* -------------------------------------------------------------------- */
147  /* Metadata. */
148  /* -------------------------------------------------------------------- */
149  char** md;
150  md = GDALGetMetadata(hDS, NULL);
151  if (md != NULL) translate_metadata(md, attr_table);
152 
153  /* -------------------------------------------------------------------- */
154  /* SRS */
155  /* -------------------------------------------------------------------- */
156  const char* pszWKT = GDALGetProjectionRef(hDS);
157  if (pszWKT != NULL && strlen(pszWKT) > 0) attach_str_attr_item(attr_table, "spatial_ref", pszWKT);
158 }
159 
167 static void build_variable_attributes(const GDALDatasetH &hDS, AttrTable *band_attr, const int iBand)
168 {
169  /* -------------------------------------------------------------------- */
170  /* Offset. */
171  /* -------------------------------------------------------------------- */
172  int bSuccess;
173  double dfValue;
174  char szValue[128];
175  GDALRasterBandH hBand = GDALGetRasterBand(hDS, iBand + 1);
176 
177  dfValue = GDALGetRasterOffset(hBand, &bSuccess);
178  if (bSuccess) {
179  snprintf(szValue, 128, "%.16g", dfValue);
180  band_attr->append_attr("add_offset", "Float64", szValue);
181  }
182 
183  /* -------------------------------------------------------------------- */
184  /* Scale */
185  /* -------------------------------------------------------------------- */
186  dfValue = GDALGetRasterScale(hBand, &bSuccess);
187  if (bSuccess) {
188  snprintf(szValue, 128, "%.16g", dfValue);
189  band_attr->append_attr("scale_factor", "Float64", szValue);
190  }
191 
192  /* -------------------------------------------------------------------- */
193  /* nodata/missing_value */
194  /* -------------------------------------------------------------------- */
195  dfValue = GDALGetRasterNoDataValue(hBand, &bSuccess);
196  if (bSuccess) {
197  snprintf(szValue, 128, "%.16g", dfValue);
198  band_attr->append_attr("missing_value", "Float64", szValue);
199  }
200 
201  /* -------------------------------------------------------------------- */
202  /* Description. */
203  /* -------------------------------------------------------------------- */
204  if (GDALGetDescription(hBand) != NULL && strlen(GDALGetDescription(hBand)) > 0) {
205  attach_str_attr_item(band_attr, "Description", GDALGetDescription(hBand));
206  }
207 
208  /* -------------------------------------------------------------------- */
209  /* PhotometricInterpretation. */
210  /* -------------------------------------------------------------------- */
211  if (GDALGetRasterColorInterpretation(hBand) != GCI_Undefined) {
212  attach_str_attr_item(band_attr, "PhotometricInterpretation",
213  GDALGetColorInterpretationName(GDALGetRasterColorInterpretation(hBand)));
214  }
215 
216  /* -------------------------------------------------------------------- */
217  /* Band Metadata. */
218  /* -------------------------------------------------------------------- */
219  char **md = GDALGetMetadata(hBand, NULL);
220  if (md != NULL) translate_metadata(md, band_attr);
221 
222  /* -------------------------------------------------------------------- */
223  /* Colormap. */
224  /* -------------------------------------------------------------------- */
225  GDALColorTableH hCT;
226 
227  hCT = GDALGetRasterColorTable(hBand);
228  if (hCT != NULL) {
229  AttrTable *ct_attr;
230  int iColor, nColorCount = GDALGetColorEntryCount(hCT);
231 
232  ct_attr = band_attr->append_container(string("Colormap"));
233 
234  for (iColor = 0; iColor < nColorCount; iColor++) {
235  GDALColorEntry sRGB;
236  AttrTable *color_attr;
237 
238  color_attr = ct_attr->append_container(string(CPLSPrintf("color_%d", iColor)));
239 
240  GDALGetColorEntryAsRGB(hCT, iColor, &sRGB);
241 
242  color_attr->append_attr("red", "byte", CPLSPrintf("%d", sRGB.c1));
243  color_attr->append_attr("green", "byte", CPLSPrintf("%d", sRGB.c2));
244  color_attr->append_attr("blue", "byte", CPLSPrintf("%d", sRGB.c3));
245  color_attr->append_attr("alpha", "byte", CPLSPrintf("%d", sRGB.c4));
246  }
247  }
248 }
249 
263 void gdal_read_dataset_attributes(DAS &das, const GDALDatasetH &hDS)
264 {
265  AttrTable *attr_table = das.add_table(string("GLOBAL"), new AttrTable);
266 
267  build_global_attributes(hDS, attr_table);
268 
269  /* ==================================================================== */
270  /* Generation info for bands. */
271  /* ==================================================================== */
272  for (int iBand = 0; iBand < GDALGetRasterCount(hDS); iBand++) {
273  char szName[128];
274  AttrTable *band_attr;
275 
276  /* -------------------------------------------------------------------- */
277  /* Create container named after the band. */
278  /* -------------------------------------------------------------------- */
279  snprintf(szName, 128, "band_%d", iBand + 1);
280  band_attr = das.add_table(string(szName), new AttrTable);
281 
282  build_variable_attributes(hDS, band_attr, iBand);
283  }
284 }
285 
295 void gdal_read_dataset_variables(DDS *dds, const GDALDatasetH &hDS, const string &filename)
296 {
297  // Load in to global attributes
298  AttrTable *global_attr = dds->get_attr_table().append_container("GLOBAL");
299  build_global_attributes(hDS, global_attr);
300 
301  /* -------------------------------------------------------------------- */
302  /* Create the basic matrix for each band. */
303  /* -------------------------------------------------------------------- */
304  BaseTypeFactory factory; // Use this to build new scalar variables
305  GDALDataType eBufType;
306 
307  for (int iBand = 0; iBand < GDALGetRasterCount(hDS); iBand++) {
308  DBG(cerr << "In dgal_dds.cc iBand" << endl);
309 
310  GDALRasterBandH hBand = GDALGetRasterBand(hDS, iBand + 1);
311 
312  ostringstream oss;
313  oss << "band_" << iBand + 1;
314 
315  eBufType = GDALGetRasterDataType(hBand);
316 
317  BaseType *bt;
318  switch (GDALGetRasterDataType(hBand)) {
319  case GDT_Byte:
320  bt = factory.NewByte(oss.str());
321  break;
322 
323  case GDT_UInt16:
324  bt = factory.NewUInt16(oss.str());
325  break;
326 
327  case GDT_Int16:
328  bt = factory.NewInt16(oss.str());
329  break;
330 
331  case GDT_UInt32:
332  bt = factory.NewUInt32(oss.str());
333  break;
334 
335  case GDT_Int32:
336  bt = factory.NewInt32(oss.str());
337  break;
338 
339  case GDT_Float32:
340  bt = factory.NewFloat32(oss.str());
341  break;
342 
343  case GDT_Float64:
344  bt = factory.NewFloat64(oss.str());
345  break;
346 
347  case GDT_CFloat32:
348  case GDT_CFloat64:
349  case GDT_CInt16:
350  case GDT_CInt32:
351  default:
352  // TODO eventually we need to preserve complex info
353  bt = factory.NewFloat64(oss.str());
354  eBufType = GDT_Float64;
355  break;
356  }
357 
358  /* -------------------------------------------------------------------- */
359  /* Create a grid to hold the raster. */
360  /* -------------------------------------------------------------------- */
361  Grid *grid;
362  grid = new GDALGrid(filename, oss.str());
363 
364  /* -------------------------------------------------------------------- */
365  /* Make into an Array for the raster data with appropriate */
366  /* dimensions. */
367  /* -------------------------------------------------------------------- */
368  Array *ar;
369  // A 'feature' of Array is that it copies the variable passed to
370  // its ctor. To get around that, pass null and use add_var_nocopy().
371  // Modified for the DAP4 response; switched to this new ctor.
372  ar = new GDALArray(oss.str(), 0, filename, eBufType, iBand + 1);
373  ar->add_var_nocopy(bt);
374  ar->append_dim(GDALGetRasterYSize(hDS), "northing");
375  ar->append_dim(GDALGetRasterXSize(hDS), "easting");
376 
377  grid->add_var_nocopy(ar, libdap::array);
378 
379  /* -------------------------------------------------------------------- */
380  /* Add the dimension map arrays. */
381  /* -------------------------------------------------------------------- */
382  //bt = new GDALFloat64( "northing" );
383  bt = factory.NewFloat64("northing");
384  ar = new GDALArray("northing", 0, filename, GDT_Float64/* eBufType */, iBand + 1);
385  ar->add_var_nocopy(bt);
386  ar->append_dim(GDALGetRasterYSize(hDS), "northing");
387 
388  grid->add_var_nocopy(ar, maps);
389 
390  //bt = new GDALFloat64( "easting" );
391  bt = factory.NewFloat64("easting");
392  ar = new GDALArray("easting", 0, filename, GDT_Float64/* eBufType */, iBand + 1);
393  ar->add_var_nocopy(bt);
394  ar->append_dim(GDALGetRasterXSize(hDS), "easting");
395 
396  grid->add_var_nocopy(ar, maps);
397 
398  DBG(cerr << "Type of grid: " << typeid(grid).name() << endl);
399 
400  // Add attributes to the Grid
401 
402  AttrTable &band_attr = grid->get_attr_table();
403  build_variable_attributes(hDS, &band_attr, iBand);
404 
405  dds->add_var_nocopy(grid);
406  }
407 }
408 
422 void gdal_read_dataset_variables(DMR *dmr, const GDALDatasetH &hDS, const string &filename)
423 {
424  // Load in global attributes. Hack, use DAP2 attributes and transform.
425  // This is easier than writing new D4Attributes code and not a heuristic
426  // routine like the variable transforms or attribute to DDS merge. New
427  // code for the D4Attributes would duplicate the AttrTable code.
428  //
429  // An extra AttrTable is needed because of oddities of its API but that
430  // comes in handy since D4Attributes::transform_to_dap4() throws away
431  // the top most AttrTable
432 
433  AttrTable *attr = new AttrTable();
434  AttrTable *global_attr = attr->append_container("GLOBAL");
435 
436  build_global_attributes(hDS, global_attr);
437 
438  dmr->root()->attributes()->transform_to_dap4(*attr);
439  delete attr; attr = 0;
440 
441  D4BaseTypeFactory factory; // Use this to build new variables
442 
443  // Make the northing and easting shared dims for this gdal dataset.
444  // For bands that have a different size than the overall Raster{X,Y}Size,
445  // assume they are lower resolution bands. For these I'm going to simply
446  // not include the shared dimensions for them. If this is a problem,
447  // we can expand the set of dimensions later.
448  //
449  // See the GDAL data model doc (http://www.gdal.org/gdal_datamodel.html)
450  // for info on why this seems correct. jhrg 6/2/17
451 
452  // Find the first band that is the same size as the whole raster dataset (i.e.,
453  // the first band that is not at a reduced resolution). In the larger loop
454  // below we only bind sdims to bands that match the overall raster size and
455  // leave bands that are at a reduce resolution w/o shared dims.
456  int sdim_band_num = 1;
457  for (; sdim_band_num <= GDALGetRasterCount(hDS); ++sdim_band_num) {
458  if (GDALGetRasterBandYSize(GDALGetRasterBand(hDS, sdim_band_num)) == GDALGetRasterYSize(hDS)
459  && GDALGetRasterBandXSize(GDALGetRasterBand(hDS, sdim_band_num)) == GDALGetRasterXSize(hDS)) {
460  break;
461  }
462  }
463 
464  // Make the two shared dims that will have the geo-referencing info
465  const string northing = "northing";
466  // Add the shared dim
467  D4Dimension *northing_dim = new D4Dimension(northing, GDALGetRasterYSize(hDS));
468  dmr->root()->dims()->add_dim_nocopy(northing_dim);
469  // use the shared dim for the map
470  Array *northing_map = new GDALArray(northing, 0, filename, GDT_Float64, sdim_band_num);
471  northing_map->add_var_nocopy(factory.NewFloat64(northing));
472  northing_map->append_dim(northing_dim);
473  // add the map
474  dmr->root()->add_var_nocopy(northing_map); // Add the map array to the DMR/D4Group
475 
476  const string easting = "easting";
477  D4Dimension *easting_dim = new D4Dimension(easting, GDALGetRasterXSize(hDS));
478  dmr->root()->dims()->add_dim_nocopy(easting_dim);
479  Array *easting_map = new GDALArray(easting, 0, filename, GDT_Float64, sdim_band_num);
480  easting_map->add_var_nocopy(factory.NewFloat64(easting));
481  easting_map->append_dim(easting_dim);
482  dmr->root()->add_var_nocopy(easting_map); // Add the map array to the DMR/D4Group
483 
484  /* -------------------------------------------------------------------- */
485  /* Create the basic matrix for each band. */
486  /* -------------------------------------------------------------------- */
487 
488  GDALDataType eBufType;
489 
490  for (int iBand = 0; iBand < GDALGetRasterCount(hDS); iBand++) {
491  DBG(cerr << "In dgal_dds.cc iBand" << endl);
492 
493  GDALRasterBandH hBand = GDALGetRasterBand(hDS, iBand + 1);
494 
495  ostringstream oss;
496  oss << "band_" << iBand + 1;
497 
498  eBufType = GDALGetRasterDataType(hBand);
499 
500  BaseType *bt;
501  switch (GDALGetRasterDataType(hBand)) {
502  case GDT_Byte:
503  bt = factory.NewByte(oss.str());
504  break;
505 
506  case GDT_UInt16:
507  bt = factory.NewUInt16(oss.str());
508  break;
509 
510  case GDT_Int16:
511  bt = factory.NewInt16(oss.str());
512  break;
513 
514  case GDT_UInt32:
515  bt = factory.NewUInt32(oss.str());
516  break;
517 
518  case GDT_Int32:
519  bt = factory.NewInt32(oss.str());
520  break;
521 
522  case GDT_Float32:
523  bt = factory.NewFloat32(oss.str());
524  break;
525 
526  case GDT_Float64:
527  bt = factory.NewFloat64(oss.str());
528  break;
529 
530  case GDT_CFloat32:
531  case GDT_CFloat64:
532  case GDT_CInt16:
533  case GDT_CInt32:
534  default:
535  // TODO eventually we need to preserve complex info
536  // Replace this with an attribute about a missing/elided variable?
537  bt = factory.NewFloat64(oss.str());
538  eBufType = GDT_Float64;
539  break;
540  }
541 
542  // Make the array for this band and then associate the dimensions/maps
543  // once they are made;
544  Array *ar = new GDALArray(oss.str(), 0, filename, eBufType, iBand + 1);
545  ar->add_var_nocopy(bt); // bt is from the above switch
546 
547  /* -------------------------------------------------------------------- */
548  /* Add the dimension map arrays. */
549  /* -------------------------------------------------------------------- */
550 
551  if (GDALGetRasterBandYSize(hBand) == GDALGetRasterYSize(hDS)
552  && GDALGetRasterBandXSize(hBand) == GDALGetRasterXSize(hDS)) {
553  // Use the shared dimensions
554  ar->append_dim(northing_dim);
555  ar->append_dim(easting_dim);
556 
557  // Bind the map to the array; FQNs for the maps (shared dims)
558  ar->maps()->add_map(new D4Map(string("/") + northing, northing_map, ar));
559  ar->maps()->add_map(new D4Map(string("/") + easting, easting_map, ar));
560  }
561  else {
562  // Don't use the shared dims
563  ar->append_dim(GDALGetRasterBandYSize(hBand), northing);
564  ar->append_dim(GDALGetRasterBandXSize(hBand), easting);
565  }
566 
567  // Add attributes, using the same hack as for the global attributes;
568  // build the DAP2 AttrTable object and then transform to DAP4.
569  AttrTable &band_attr = ar->get_attr_table();
570  build_variable_attributes(hDS, &band_attr, iBand);
571  ar->attributes()->transform_to_dap4(band_attr);
572 
573  dmr->root()->add_var_nocopy(ar); // Add the array to the DMR
574  }
575 }
576 
594 void read_data_array(GDALArray *array, const GDALRasterBandH &hBand)
595 {
596  /* -------------------------------------------------------------------- */
597  /* Collect the x and y sampling values from the constraint. */
598  /* -------------------------------------------------------------------- */
599  Array::Dim_iter p = array->dim_begin();
600  int start = array->dimension_start(p, true);
601  int stride = array->dimension_stride(p, true);
602  int stop = array->dimension_stop(p, true);
603 
604  // Test for the case where a dimension has not been subset. jhrg 2/18/16
605  if (array->dimension_size(p, true) == 0) { //default rows
606  start = 0;
607  stride = 1;
608  stop = GDALGetRasterBandYSize(hBand) - 1;
609  }
610 
611  p++;
612  int start_2 = array->dimension_start(p, true);
613  int stride_2 = array->dimension_stride(p, true);
614  int stop_2 = array->dimension_stop(p, true);
615 
616  if (array->dimension_size(p, true) == 0) { //default columns
617  start_2 = 0;
618  stride_2 = 1;
619  stop_2 = GDALGetRasterBandXSize(hBand) - 1;
620  }
621 
622  /* -------------------------------------------------------------------- */
623  /* Build a window and buf size from this. */
624  /* -------------------------------------------------------------------- */
625  int nWinXOff, nWinYOff, nWinXSize, nWinYSize, nBufXSize, nBufYSize;
626 
627  nWinXOff = start_2;
628  nWinYOff = start;
629  nWinXSize = stop_2 + 1 - start_2;
630  nWinYSize = stop + 1 - start;
631 
632  nBufXSize = (stop_2 - start_2) / stride_2 + 1;
633  nBufYSize = (stop - start) / stride + 1;
634 
635  /* -------------------------------------------------------------------- */
636  /* Allocate buffer. */
637  /* -------------------------------------------------------------------- */
638  int nPixelSize = GDALGetDataTypeSize(array->get_gdal_buf_type()) / 8;
639  vector<char> pData(nBufXSize * nBufYSize * nPixelSize);
640 
641  /* -------------------------------------------------------------------- */
642  /* Read request into buffer. */
643  /* -------------------------------------------------------------------- */
644  CPLErr eErr = GDALRasterIO(hBand, GF_Read, nWinXOff, nWinYOff, nWinXSize, nWinYSize, &pData[0], nBufXSize,
645  nBufYSize, array->get_gdal_buf_type(), 0, 0);
646  if (eErr != CE_None) throw Error("Error reading: " + array->name());
647 
648  array->val2buf(&pData[0]);
649 }
650 
669 void read_map_array(Array *map, const GDALRasterBandH &hBand, const GDALDatasetH &hDS)
670 {
671  Array::Dim_iter p = map->dim_begin();
672  int start = map->dimension_start(p, true);
673  int stride = map->dimension_stride(p, true);
674  int stop = map->dimension_stop(p, true);
675 
676  if (start + stop + stride == 0) { //default rows
677  start = 0;
678  stride = 1;
679  if (map->name() == "northing")
680  stop = GDALGetRasterBandYSize(hBand) - 1;
681  else if (map->name() == "easting")
682  stop = GDALGetRasterBandXSize(hBand) - 1;
683  else
684  throw Error("Expected a map named 'northing' or 'easting' but got: " + map->name());
685  }
686 
687  int nBufSize = (stop - start) / stride + 1;
688 
689  /* -------------------------------------------------------------------- */
690  /* Read or default the geotransform used to generate the */
691  /* georeferencing maps. */
692  /* -------------------------------------------------------------------- */
693 
694  double adfGeoTransform[6];
695 
696  if (GDALGetGeoTransform(hDS, adfGeoTransform) != CE_None) {
697  adfGeoTransform[0] = 0.0;
698  adfGeoTransform[1] = 1.0;
699  adfGeoTransform[2] = 0.0;
700  adfGeoTransform[3] = 0.0;
701  adfGeoTransform[4] = 0.0;
702  adfGeoTransform[5] = 1.0;
703  }
704 
705  /* -------------------------------------------------------------------- */
706  /* Set the map array. */
707  /* -------------------------------------------------------------------- */
708  vector<double> padfMap(nBufSize);
709 
710  if (map->name() == "northing") {
711  for (int i = 0, iLine = start; iLine <= stop; iLine += stride) {
712  padfMap[i++] = adfGeoTransform[3] + adfGeoTransform[5] * iLine;
713  }
714  }
715  else if (map->name() == "easting") {
716  for (int i = 0, iPixel = start; iPixel <= stop; iPixel += stride) {
717  padfMap[i++] = adfGeoTransform[0] + iPixel * adfGeoTransform[1];
718  }
719  }
720  else
721  throw Error("Expected a map named 'northing' or 'easting' but got: " + map->name());
722 
723  map->val2buf((void *) &padfMap[0]);
724 }
725 
726 // $Log: gdal_das.cc,v $
727 // Revision 1.1 2004/10/19 20:38:28 warmerda
728 // New
729 //
730 // Revision 1.2 2004/10/15 18:06:45 warmerda
731 // Strip the extension off the filename.
732 //
733 // Revision 1.1 2004/10/04 14:29:29 warmerda
734 // New
735 //
736 
libdap
Definition: BESDapFunctionResponseCache.h:35
GDALGrid
Definition: GDALTypes.h:59
GDALArray
Definition: GDALTypes.h:34
Error