bes Updated for version 3.20.13
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 "GDALTypes.h"
54
55using 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
64static 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
82static 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
107static 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
171static 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
267void 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
300void 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
430void 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
602void 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.data(), 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.data());
657}
658
677void 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.data());
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