29#include <cpl_string.h>
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>
41#include <libdap/DDS.h>
42#include <libdap/DAS.h>
43#include <libdap/BaseTypeFactory.h>
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>
64static void attach_str_attr_item(AttrTable *parent_table,
const char *pszKey,
const char *pszValue)
67 char *pszEscapedText = CPLEscapeString(pszValue, -1,
68 CPLES_BackslashQuotable);
70 parent_table->append_attr(pszKey,
"String", pszEscapedText );
72 CPLFree(pszEscapedText);
82static void translate_metadata(
char **md, AttrTable *parent_table)
87 md_table = parent_table->append_container(
string(
"Metadata"));
89 for (i = 0; md != NULL && md[i] != NULL; i++) {
93 pszValue = CPLParseNameValue(md[i], &pszKey);
95 attach_str_attr_item(md_table, pszKey, pszValue);
107static void build_global_attributes(
const GDALDatasetH& hDS, AttrTable* attr_table)
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)) {
117 char szGeoTransform[200];
118 double dfMaxX, dfMinX, dfMaxY, dfMinY;
119 int nXSize = GDALGetRasterXSize(hDS);
120 int nYSize = GDALGetRasterYSize(hDS);
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));
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));
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));
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));
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));
140 attr_table->append_attr(
"Westernmost_Northing",
"Float64", CPLSPrintf(
"%.16g", dfMinX));
142 attr_table->append_attr(
"Westernmost_Easting",
"Float64", CPLSPrintf(
"%.16g", dfMinX));
144 snprintf(szGeoTransform, 200,
"%.16g %.16g %.16g %.16g %.16g %.16g", adfGeoTransform[0], adfGeoTransform[1],
145 adfGeoTransform[2], adfGeoTransform[3], adfGeoTransform[4], adfGeoTransform[5]);
147 attach_str_attr_item(attr_table,
"GeoTransform", szGeoTransform);
154 md = GDALGetMetadata(hDS, NULL);
155 if (md != NULL) translate_metadata(md, attr_table);
160 const char* pszWKT = GDALGetProjectionRef(hDS);
161 if (pszWKT != NULL && strlen(pszWKT) > 0) attach_str_attr_item(attr_table,
"spatial_ref", pszWKT);
171static void build_variable_attributes(
const GDALDatasetH &hDS, AttrTable *band_attr,
const int iBand)
179 GDALRasterBandH hBand = GDALGetRasterBand(hDS, iBand + 1);
181 dfValue = GDALGetRasterOffset(hBand, &bSuccess);
183 snprintf(szValue, 128,
"%.16g", dfValue);
184 band_attr->append_attr(
"add_offset",
"Float64", szValue);
190 dfValue = GDALGetRasterScale(hBand, &bSuccess);
192 snprintf(szValue, 128,
"%.16g", dfValue);
193 band_attr->append_attr(
"scale_factor",
"Float64", szValue);
199 dfValue = GDALGetRasterNoDataValue(hBand, &bSuccess);
201 snprintf(szValue, 128,
"%.16g", dfValue);
202 band_attr->append_attr(
"missing_value",
"Float64", szValue);
208 if (GDALGetDescription(hBand) != NULL && strlen(GDALGetDescription(hBand)) > 0) {
209 attach_str_attr_item(band_attr,
"Description", GDALGetDescription(hBand));
215 if (GDALGetRasterColorInterpretation(hBand) != GCI_Undefined) {
216 attach_str_attr_item(band_attr,
"PhotometricInterpretation",
217 GDALGetColorInterpretationName(GDALGetRasterColorInterpretation(hBand)));
223 char **md = GDALGetMetadata(hBand, NULL);
224 if (md != NULL) translate_metadata(md, band_attr);
231 hCT = GDALGetRasterColorTable(hBand);
234 int iColor, nColorCount = GDALGetColorEntryCount(hCT);
236 ct_attr = band_attr->append_container(
string(
"Colormap"));
238 for (iColor = 0; iColor < nColorCount; iColor++) {
240 AttrTable *color_attr;
242 color_attr = ct_attr->append_container(
string(CPLSPrintf(
"color_%d", iColor)));
244 GDALGetColorEntryAsRGB(hCT, iColor, &sRGB);
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));
267void gdal_read_dataset_attributes(DAS &das,
const GDALDatasetH &hDS)
269 AttrTable *attr_table = das.add_table(
string(
"GLOBAL"),
new AttrTable);
271 build_global_attributes(hDS, attr_table);
276 for (
int iBand = 0; iBand < GDALGetRasterCount(hDS); iBand++) {
278 AttrTable *band_attr;
283 snprintf(szName, 128,
"band_%d", iBand + 1);
284 band_attr = das.add_table(
string(szName),
new AttrTable);
286 build_variable_attributes(hDS, band_attr, iBand);
300void gdal_read_dataset_variables(DDS *dds,
const GDALDatasetH &hDS,
const string &filename,
bool include_attrs)
303 if(
true == include_attrs) {
304 AttrTable *global_attr = dds->get_attr_table().append_container(
"GLOBAL");
305 build_global_attributes(hDS, global_attr);
312 GDALDataType eBufType;
314 for (
int iBand = 0; iBand < GDALGetRasterCount(hDS); iBand++) {
315 DBG(cerr <<
"In dgal_dds.cc iBand" << endl);
317 GDALRasterBandH hBand = GDALGetRasterBand(hDS, iBand + 1);
320 oss <<
"band_" << iBand + 1;
322 eBufType = GDALGetRasterDataType(hBand);
325 switch (GDALGetRasterDataType(hBand)) {
327 bt = factory.NewByte(oss.str());
331 bt = factory.NewUInt16(oss.str());
335 bt = factory.NewInt16(oss.str());
339 bt = factory.NewUInt32(oss.str());
343 bt = factory.NewInt32(oss.str());
347 bt = factory.NewFloat32(oss.str());
351 bt = factory.NewFloat64(oss.str());
360 bt = factory.NewFloat64(oss.str());
361 eBufType = GDT_Float64;
369 grid =
new GDALGrid(filename, oss.str());
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");
384 grid->add_var_nocopy(ar, libdap::array);
390 bt = factory.NewFloat64(
"northing");
391 ar =
new GDALArray(
"northing", 0, filename, GDT_Float64, iBand + 1);
392 ar->add_var_nocopy(bt);
393 ar->append_dim(GDALGetRasterYSize(hDS),
"northing");
395 grid->add_var_nocopy(ar, maps);
398 bt = factory.NewFloat64(
"easting");
399 ar =
new GDALArray(
"easting", 0, filename, GDT_Float64, iBand + 1);
400 ar->add_var_nocopy(bt);
401 ar->append_dim(GDALGetRasterXSize(hDS),
"easting");
403 grid->add_var_nocopy(ar, maps);
405 DBG(cerr <<
"Type of grid: " <<
typeid(grid).name() << endl);
408 if(
true == include_attrs) {
409 AttrTable &band_attr = grid->get_attr_table();
410 build_variable_attributes(hDS, &band_attr, iBand);
413 dds->add_var_nocopy(grid);
430void gdal_read_dataset_variables(DMR *dmr,
const GDALDatasetH &hDS,
const string &filename)
441 AttrTable *attr =
new AttrTable();
442 AttrTable *global_attr = attr->append_container(
"GLOBAL");
444 build_global_attributes(hDS, global_attr);
446 dmr->root()->attributes()->transform_to_dap4(*attr);
447 delete attr; attr = 0;
449 D4BaseTypeFactory factory;
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)) {
473 const string northing =
"northing";
475 D4Dimension *northing_dim =
new D4Dimension(northing, GDALGetRasterYSize(hDS));
476 dmr->root()->dims()->add_dim_nocopy(northing_dim);
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);
482 dmr->root()->add_var_nocopy(northing_map);
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);
496 GDALDataType eBufType;
498 for (
int iBand = 0; iBand < GDALGetRasterCount(hDS); iBand++) {
499 DBG(cerr <<
"In dgal_dds.cc iBand" << endl);
501 GDALRasterBandH hBand = GDALGetRasterBand(hDS, iBand + 1);
504 oss <<
"band_" << iBand + 1;
506 eBufType = GDALGetRasterDataType(hBand);
509 switch (GDALGetRasterDataType(hBand)) {
511 bt = factory.NewByte(oss.str());
515 bt = factory.NewUInt16(oss.str());
519 bt = factory.NewInt16(oss.str());
523 bt = factory.NewUInt32(oss.str());
527 bt = factory.NewInt32(oss.str());
531 bt = factory.NewFloat32(oss.str());
535 bt = factory.NewFloat64(oss.str());
545 bt = factory.NewFloat64(oss.str());
546 eBufType = GDT_Float64;
552 Array *ar =
new GDALArray(oss.str(), 0, filename, eBufType, iBand + 1);
553 ar->add_var_nocopy(bt);
559 if (GDALGetRasterBandYSize(hBand) == GDALGetRasterYSize(hDS)
560 && GDALGetRasterBandXSize(hBand) == GDALGetRasterXSize(hDS)) {
562 ar->append_dim(northing_dim);
563 ar->append_dim(easting_dim);
566 ar->maps()->add_map(
new D4Map(
string(
"/") + northing, northing_map, ar));
567 ar->maps()->add_map(
new D4Map(
string(
"/") + easting, easting_map, ar));
571 ar->append_dim(GDALGetRasterBandYSize(hBand), northing);
572 ar->append_dim(GDALGetRasterBandXSize(hBand), easting);
577 AttrTable &band_attr = ar->get_attr_table();
578 build_variable_attributes(hDS, &band_attr, iBand);
579 ar->attributes()->transform_to_dap4(band_attr);
581 dmr->root()->add_var_nocopy(ar);
602void read_data_array(
GDALArray *array,
const GDALRasterBandH &hBand)
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);
613 if (array->dimension_size(p,
true) == 0) {
616 stop = GDALGetRasterBandYSize(hBand) - 1;
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);
624 if (array->dimension_size(p,
true) == 0) {
627 stop_2 = GDALGetRasterBandXSize(hBand) - 1;
633 int nWinXOff, nWinYOff, nWinXSize, nWinYSize, nBufXSize, nBufYSize;
637 nWinXSize = stop_2 + 1 - start_2;
638 nWinYSize = stop + 1 - start;
640 nBufXSize = (stop_2 - start_2) / stride_2 + 1;
641 nBufYSize = (stop - start) / stride + 1;
646 int nPixelSize = GDALGetDataTypeSize(array->get_gdal_buf_type()) / 8;
647 vector<char> pData(nBufXSize * nBufYSize * nPixelSize);
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());
656 array->val2buf(pData.data());
677void read_map_array(Array *map,
const GDALRasterBandH &hBand,
const GDALDatasetH &hDS)
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);
684 if (start + stop + stride == 0) {
687 if (map->name() ==
"northing")
688 stop = GDALGetRasterBandYSize(hBand) - 1;
689 else if (map->name() ==
"easting")
690 stop = GDALGetRasterBandXSize(hBand) - 1;
692 throw Error(
"Expected a map named 'northing' or 'easting' but got: " + map->name());
695 int nBufSize = (stop - start) / stride + 1;
702 double adfGeoTransform[6];
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;
716 vector<double> padfMap(nBufSize);
718 if (map->name() ==
"northing") {
719 for (
int i = 0, iLine = start; iLine <= stop; iLine += stride) {
720 padfMap[i++] = adfGeoTransform[3] + adfGeoTransform[5] * iLine;
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];
729 throw Error(
"Expected a map named 'northing' or 'easting' but got: " + map->name());
731 map->val2buf((
void *) padfMap.data());