13 #include "SpatialIndex.h"
14 #include "SpatialInterface.h"
16 #include <D4Connect.h>
23 #include <unordered_map>
30 static bool verbose =
false;
31 #define VERBOSE(x) do { if (verbose) x; } while(false)
51 void findLatLon(std::string dataUrl,
const float64 level,
52 const float64 buildlevel, vector<int> &xArray, vector<int> &yArray, vector<float> &latArray, vector<float> &lonArray, vector<uint64> &stareArray) {
54 htmInterface htm(level, buildlevel);
55 const SpatialIndex &index = htm.index();
57 auto_ptr < libdap::Connect > url(
new libdap::Connect(dataUrl));
59 std::vector<float> lat;
60 std::vector<float> lon;
63 libdap::BaseTypeFactory factory;
64 libdap::DataDDS dds(&factory);
66 VERBOSE(cerr <<
"\n\n\tRequesting data from " << dataUrl << endl);
69 url->request_data(dds,
"Latitude,Longitude");
72 libdap::Array *urlLatArray = dynamic_cast<libdap::Array *>(dds.var(
74 libdap::Array *urlLonArray = dynamic_cast<libdap::Array *>(dds.var(
78 if (urlLatArray == 0 || urlLonArray == 0) {
79 throw libdap::Error(
"Expected both lat and lon arrays");
82 unsigned int dims = urlLatArray->dimensions();
85 throw libdap::Error(
"Incorrect latitude dimensions");
88 if (dims != urlLonArray->dimensions()) {
89 throw libdap::Error(
"Incorrect longitude dimensions");
92 int size_y = urlLatArray->dimension_size(urlLatArray->dim_begin());
93 int size_x = urlLatArray->dimension_size(urlLatArray->dim_begin() + 1);
95 if (size_y != urlLonArray->dimension_size(urlLonArray->dim_begin())
97 != urlLonArray->dimension_size(
98 urlLonArray->dim_begin() + 1)) {
100 "The size of the latitude and longitude arrays are not the same");
104 lat.resize(urlLatArray->length());
105 urlLatArray->value(&lat[0]);
106 lon.resize(urlLonArray->length());
107 urlLonArray->value(&lon[0]);
109 VERBOSE(cerr <<
"\tsize of lat array: " << lat.size() << endl);
110 VERBOSE(cerr <<
"\tsize of lon array: " << lon.size() << endl);
112 coords indexVals = coords();
113 std::unordered_map<float, struct coords> indexMap;
121 keyVals.resize(size_x * size_y);
124 xArray.resize(lat.size());
125 yArray.resize(lat.size());
126 latArray.resize(lat.size());
127 lonArray.resize(lat.size());
128 stareArray.resize(lat.size());
132 vector<float>::iterator j_begin = lon.begin();
136 VERBOSE (cerr <<
"\nCalculating the STARE indices" << endl);
140 for (vector<float>::iterator i = lat.begin(), e = lat.end(), j =
141 lon.begin(); i != e; ++i, ++j) {
143 int offset = j - j_begin;
144 indexVals.x = offset / (size_x - 1);
145 indexVals.y = offset % (size_x - 1);
148 indexVals.stareIndex = index.idByLatLon(*i, *j);
151 indexMap[indexVals.stareIndex] = indexVals;
154 xArray[arrayLoc] = indexVals.x;
155 yArray[arrayLoc] = indexVals.y;
156 latArray[arrayLoc] = *i;
157 lonArray[arrayLoc] = *j;
158 stareArray[arrayLoc] = indexVals.stareIndex;
166 keyVals[arrayLoc].x = indexVals.x;
167 keyVals[arrayLoc].y = indexVals.y;
168 keyVals[arrayLoc].lat = *i;
169 keyVals[arrayLoc].lon = *j;
170 keyVals[arrayLoc].stareIndex = indexVals.stareIndex;
176 }
catch (libdap::Error &e) {
177 cerr <<
"ERROR: " << e.get_error_message() << endl;
185 void writeHDF5(
const string &filename,
const vector<int> &xArray,
const vector<int> &yArray,
const vector<float> &latArray,
const vector<float> &lonArray,
const vector<uint64> &stareArray) {
186 hid_t file, datasetX, datasetY, datasetLat, datasetLon, datasetStare;
191 hsize_t arrayLength[1] = { xArray.size() };
196 hid_t fapl = H5Pcreate (H5P_FILE_ACCESS);
197 #ifdef H5F_LIBVER_V18
198 H5Pset_libver_bounds (fapl, H5F_LIBVER_EARLIEST, H5F_LIBVER_V18);
200 H5Pset_libver_bounds (fapl, H5F_LIBVER_EARLIEST, H5F_LIBVER_LATEST);
204 file = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, fapl);
208 dataspace = H5Screate_simple(1 , arrayLength, NULL);
220 VERBOSE(cerr <<
"\nCreating datatypes: x, y, lat, lon, stareIndex -> ");
222 dataTypes = H5Tcreate (H5T_COMPOUND,
sizeof(coords));
223 H5Tinsert(dataTypes,
"x", HOFFSET(coords, x), H5T_NATIVE_INT);
224 H5Tinsert(dataTypes,
"y", HOFFSET(coords, y), H5T_NATIVE_INT);
225 H5Tinsert(dataTypes,
"lat", HOFFSET(coords, lat), H5T_NATIVE_FLOAT);
226 H5Tinsert(dataTypes,
"lon", HOFFSET(coords, lon), H5T_NATIVE_FLOAT);
227 H5Tinsert(dataTypes,
"stareIndex", HOFFSET(coords, stareIndex), H5T_NATIVE_INT);
232 const char *datasetName =
"StareData";
233 VERBOSE(cerr <<
"Creating dataset: " << datasetName <<
" -> ");
234 dataset = H5Dcreate2(file, datasetName, dataTypes, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
236 VERBOSE(cerr <<
"Writing data to dataset" << endl);
237 H5Dwrite(dataset, dataTypes, H5S_ALL, H5S_ALL, H5P_DEFAULT, &keyVals[0]);
243 const char *datasetName =
"X";
244 VERBOSE(cerr <<
"Creating dataset: " << datasetName <<
" -> ");
245 datasetX = H5Dcreate2(file, datasetName, H5T_NATIVE_INT, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
248 VERBOSE(cerr <<
"Creating dataset: " << datasetName <<
" -> ");
249 datasetY = H5Dcreate2(file, datasetName, H5T_NATIVE_INT, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
251 datasetName =
"Latitude";
252 VERBOSE(cerr <<
"Creating dataset: " << datasetName <<
" -> ");
253 datasetLat = H5Dcreate2(file, datasetName, H5T_NATIVE_FLOAT, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
255 datasetName =
"Longitude";
256 VERBOSE(cerr <<
"Creating dataset: " << datasetName <<
" -> ");
257 datasetLon = H5Dcreate2(file, datasetName, H5T_NATIVE_FLOAT, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
259 datasetName =
"Stare Index";
260 VERBOSE(cerr <<
"Creating dataset: " << datasetName <<
" -> ");
261 datasetStare = H5Dcreate2(file, datasetName, H5T_NATIVE_INT64, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
266 VERBOSE(cerr <<
"Writing data to dataset" << endl);
267 H5Dwrite(datasetX, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &xArray[0]);
268 H5Dwrite(datasetY, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &yArray[0]);
269 H5Dwrite(datasetLat, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &latArray[0]);
270 H5Dwrite(datasetLon, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &lonArray[0]);
271 H5Dwrite(datasetStare, H5T_NATIVE_INT64, H5S_ALL, H5S_ALL, H5P_DEFAULT, &stareArray[0]);
279 VERBOSE(cerr <<
"Data written to file: " << filename << endl);
289 int main(
int argc,
char *argv[]) {
293 while ((c = getopt(argc, argv,
"hvo:")) != -1) {
296 cerr <<
"\nbuild_sidecar [options] <filename> level buildlevel\n\n";
297 cerr <<
"-o output file: \tOutput the STARE data to the given output file\n\n";
298 cerr <<
"-v verbose\n" << endl;
311 string dataUrl = argv[argc-3];
312 float level = atof(argv[argc-2]);
313 float build_level = atof(argv[argc-1]);
324 vector<float> latVals;
325 vector<float> lonVals;
326 vector<uint64> stareVals;
328 findLatLon(dataUrl, level, build_level, xVals, yVals, latVals, lonVals, stareVals);
334 vector<coords> keyVals;
335 findLatLon(dataUrl, level, build_level, keyVals);
338 if (newName.empty()) {
342 size_t granulePos = dataUrl.find_last_of(
"/");
343 string granuleName = dataUrl.substr(granulePos + 1);
344 size_t findDot = granuleName.find_last_of(
".");
345 newName = granuleName.substr(0, findDot) +
"_sidecar.h5";
348 writeHDF5(newName, xVals, yVals, latVals, lonVals, stareVals);
350 catch(libdap::Error &e) {
351 cerr <<
"Error: " << e.get_error_message() << endl;