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