29 #include <cpl_string.h>
43 #include <BaseTypeFactory.h>
47 #include <D4Dimensions.h>
49 #include <D4Attributes.h>
50 #include <D4BaseTypeFactory.h>
53 #include "GDALTypes.h"
64 static 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);
82 static 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);
107 static 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));
138 attr_table->append_attr(
"Westernmost_Northing",
"Float64", CPLSPrintf(
"%.16g", dfMinX));
140 snprintf(szGeoTransform, 200,
"%.16g %.16g %.16g %.16g %.16g %.16g", adfGeoTransform[0], adfGeoTransform[1],
141 adfGeoTransform[2], adfGeoTransform[3], adfGeoTransform[4], adfGeoTransform[5]);
143 attach_str_attr_item(attr_table,
"GeoTransform", szGeoTransform);
150 md = GDALGetMetadata(hDS, NULL);
151 if (md != NULL) translate_metadata(md, attr_table);
156 const char* pszWKT = GDALGetProjectionRef(hDS);
157 if (pszWKT != NULL && strlen(pszWKT) > 0) attach_str_attr_item(attr_table,
"spatial_ref", pszWKT);
167 static void build_variable_attributes(
const GDALDatasetH &hDS, AttrTable *band_attr,
const int iBand)
175 GDALRasterBandH hBand = GDALGetRasterBand(hDS, iBand + 1);
177 dfValue = GDALGetRasterOffset(hBand, &bSuccess);
179 snprintf(szValue, 128,
"%.16g", dfValue);
180 band_attr->append_attr(
"add_offset",
"Float64", szValue);
186 dfValue = GDALGetRasterScale(hBand, &bSuccess);
188 snprintf(szValue, 128,
"%.16g", dfValue);
189 band_attr->append_attr(
"scale_factor",
"Float64", szValue);
195 dfValue = GDALGetRasterNoDataValue(hBand, &bSuccess);
197 snprintf(szValue, 128,
"%.16g", dfValue);
198 band_attr->append_attr(
"missing_value",
"Float64", szValue);
204 if (GDALGetDescription(hBand) != NULL && strlen(GDALGetDescription(hBand)) > 0) {
205 attach_str_attr_item(band_attr,
"Description", GDALGetDescription(hBand));
211 if (GDALGetRasterColorInterpretation(hBand) != GCI_Undefined) {
212 attach_str_attr_item(band_attr,
"PhotometricInterpretation",
213 GDALGetColorInterpretationName(GDALGetRasterColorInterpretation(hBand)));
219 char **md = GDALGetMetadata(hBand, NULL);
220 if (md != NULL) translate_metadata(md, band_attr);
227 hCT = GDALGetRasterColorTable(hBand);
230 int iColor, nColorCount = GDALGetColorEntryCount(hCT);
232 ct_attr = band_attr->append_container(
string(
"Colormap"));
234 for (iColor = 0; iColor < nColorCount; iColor++) {
236 AttrTable *color_attr;
238 color_attr = ct_attr->append_container(
string(CPLSPrintf(
"color_%d", iColor)));
240 GDALGetColorEntryAsRGB(hCT, iColor, &sRGB);
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));
263 void gdal_read_dataset_attributes(DAS &das,
const GDALDatasetH &hDS)
265 AttrTable *attr_table = das.add_table(
string(
"GLOBAL"),
new AttrTable);
267 build_global_attributes(hDS, attr_table);
272 for (
int iBand = 0; iBand < GDALGetRasterCount(hDS); iBand++) {
274 AttrTable *band_attr;
279 snprintf(szName, 128,
"band_%d", iBand + 1);
280 band_attr = das.add_table(
string(szName),
new AttrTable);
282 build_variable_attributes(hDS, band_attr, iBand);
295 void gdal_read_dataset_variables(DDS *dds,
const GDALDatasetH &hDS,
const string &filename)
298 AttrTable *global_attr = dds->get_attr_table().append_container(
"GLOBAL");
299 build_global_attributes(hDS, global_attr);
304 BaseTypeFactory factory;
305 GDALDataType eBufType;
307 for (
int iBand = 0; iBand < GDALGetRasterCount(hDS); iBand++) {
308 DBG(cerr <<
"In dgal_dds.cc iBand" << endl);
310 GDALRasterBandH hBand = GDALGetRasterBand(hDS, iBand + 1);
313 oss <<
"band_" << iBand + 1;
315 eBufType = GDALGetRasterDataType(hBand);
318 switch (GDALGetRasterDataType(hBand)) {
320 bt = factory.NewByte(oss.str());
324 bt = factory.NewUInt16(oss.str());
328 bt = factory.NewInt16(oss.str());
332 bt = factory.NewUInt32(oss.str());
336 bt = factory.NewInt32(oss.str());
340 bt = factory.NewFloat32(oss.str());
344 bt = factory.NewFloat64(oss.str());
353 bt = factory.NewFloat64(oss.str());
354 eBufType = GDT_Float64;
362 grid =
new GDALGrid(filename, oss.str());
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");
377 grid->add_var_nocopy(ar, libdap::array);
383 bt = factory.NewFloat64(
"northing");
384 ar =
new GDALArray(
"northing", 0, filename, GDT_Float64, iBand + 1);
385 ar->add_var_nocopy(bt);
386 ar->append_dim(GDALGetRasterYSize(hDS),
"northing");
388 grid->add_var_nocopy(ar, maps);
391 bt = factory.NewFloat64(
"easting");
392 ar =
new GDALArray(
"easting", 0, filename, GDT_Float64, iBand + 1);
393 ar->add_var_nocopy(bt);
394 ar->append_dim(GDALGetRasterXSize(hDS),
"easting");
396 grid->add_var_nocopy(ar, maps);
398 DBG(cerr <<
"Type of grid: " <<
typeid(grid).name() << endl);
402 AttrTable &band_attr = grid->get_attr_table();
403 build_variable_attributes(hDS, &band_attr, iBand);
405 dds->add_var_nocopy(grid);
422 void gdal_read_dataset_variables(DMR *dmr,
const GDALDatasetH &hDS,
const string &filename)
433 AttrTable *attr =
new AttrTable();
434 AttrTable *global_attr = attr->append_container(
"GLOBAL");
436 build_global_attributes(hDS, global_attr);
438 dmr->root()->attributes()->transform_to_dap4(*attr);
439 delete attr; attr = 0;
441 D4BaseTypeFactory factory;
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)) {
465 const string northing =
"northing";
467 D4Dimension *northing_dim =
new D4Dimension(northing, GDALGetRasterYSize(hDS));
468 dmr->root()->dims()->add_dim_nocopy(northing_dim);
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);
474 dmr->root()->add_var_nocopy(northing_map);
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);
488 GDALDataType eBufType;
490 for (
int iBand = 0; iBand < GDALGetRasterCount(hDS); iBand++) {
491 DBG(cerr <<
"In dgal_dds.cc iBand" << endl);
493 GDALRasterBandH hBand = GDALGetRasterBand(hDS, iBand + 1);
496 oss <<
"band_" << iBand + 1;
498 eBufType = GDALGetRasterDataType(hBand);
501 switch (GDALGetRasterDataType(hBand)) {
503 bt = factory.NewByte(oss.str());
507 bt = factory.NewUInt16(oss.str());
511 bt = factory.NewInt16(oss.str());
515 bt = factory.NewUInt32(oss.str());
519 bt = factory.NewInt32(oss.str());
523 bt = factory.NewFloat32(oss.str());
527 bt = factory.NewFloat64(oss.str());
537 bt = factory.NewFloat64(oss.str());
538 eBufType = GDT_Float64;
544 Array *ar =
new GDALArray(oss.str(), 0, filename, eBufType, iBand + 1);
545 ar->add_var_nocopy(bt);
551 if (GDALGetRasterBandYSize(hBand) == GDALGetRasterYSize(hDS)
552 && GDALGetRasterBandXSize(hBand) == GDALGetRasterXSize(hDS)) {
554 ar->append_dim(northing_dim);
555 ar->append_dim(easting_dim);
558 ar->maps()->add_map(
new D4Map(
string(
"/") + northing, northing_map, ar));
559 ar->maps()->add_map(
new D4Map(
string(
"/") + easting, easting_map, ar));
563 ar->append_dim(GDALGetRasterBandYSize(hBand), northing);
564 ar->append_dim(GDALGetRasterBandXSize(hBand), easting);
569 AttrTable &band_attr = ar->get_attr_table();
570 build_variable_attributes(hDS, &band_attr, iBand);
571 ar->attributes()->transform_to_dap4(band_attr);
573 dmr->root()->add_var_nocopy(ar);
594 void read_data_array(
GDALArray *array,
const GDALRasterBandH &hBand)
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);
605 if (array->dimension_size(p,
true) == 0) {
608 stop = GDALGetRasterBandYSize(hBand) - 1;
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);
616 if (array->dimension_size(p,
true) == 0) {
619 stop_2 = GDALGetRasterBandXSize(hBand) - 1;
625 int nWinXOff, nWinYOff, nWinXSize, nWinYSize, nBufXSize, nBufYSize;
629 nWinXSize = stop_2 + 1 - start_2;
630 nWinYSize = stop + 1 - start;
632 nBufXSize = (stop_2 - start_2) / stride_2 + 1;
633 nBufYSize = (stop - start) / stride + 1;
638 int nPixelSize = GDALGetDataTypeSize(array->get_gdal_buf_type()) / 8;
639 vector<char> pData(nBufXSize * nBufYSize * nPixelSize);
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());
648 array->val2buf(&pData[0]);
669 void read_map_array(Array *map,
const GDALRasterBandH &hBand,
const GDALDatasetH &hDS)
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);
676 if (start + stop + stride == 0) {
679 if (map->name() ==
"northing")
680 stop = GDALGetRasterBandYSize(hBand) - 1;
681 else if (map->name() ==
"easting")
682 stop = GDALGetRasterBandXSize(hBand) - 1;
684 throw Error(
"Expected a map named 'northing' or 'easting' but got: " + map->name());
687 int nBufSize = (stop - start) / stride + 1;
694 double adfGeoTransform[6];
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;
708 vector<double> padfMap(nBufSize);
710 if (map->name() ==
"northing") {
711 for (
int i = 0, iLine = start; iLine <= stop; iLine += stride) {
712 padfMap[i++] = adfGeoTransform[3] + adfGeoTransform[5] * iLine;
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];
721 throw Error(
"Expected a map named 'northing' or 'easting' but got: " + map->name());
723 map->val2buf((
void *) &padfMap[0]);