15 #include <D4Connect.h>
22 #include <unordered_map>
29 static bool verbose =
false;
30 #define VERBOSE(x) do { if (verbose) x; } while(false)
49 vector<coords> readUrl(
string dataUrl,
string latName,
string lonName) {
50 auto_ptr < libdap::Connect > url(
new libdap::Connect(dataUrl));
52 string latlonName = latName +
"," + lonName;
54 std::vector<float> lat;
55 std::vector<float> lon;
57 vector<coords> indexArray;
60 libdap::BaseTypeFactory factory;
61 libdap::DataDDS dds(&factory);
63 VERBOSE(cerr <<
"\n\n\tRequesting data from " << dataUrl << endl);
66 url->request_data(dds, latlonName);
69 libdap::Array *urlLatArray =
dynamic_cast<libdap::Array *
>(dds.var(
71 libdap::Array *urlLonArray =
dynamic_cast<libdap::Array *
>(dds.var(
75 if (urlLatArray == 0 || urlLonArray == 0) {
76 throw libdap::Error(
"Expected both lat and lon arrays");
79 unsigned int dims = urlLatArray->dimensions();
82 throw libdap::Error(
"Incorrect latitude dimensions");
85 if (dims != urlLonArray->dimensions()) {
86 throw libdap::Error(
"Incorrect longitude dimensions");
89 int size_y = urlLatArray->dimension_size(urlLatArray->dim_begin());
90 int size_x = urlLatArray->dimension_size(urlLatArray->dim_begin() + 1);
92 if (size_y != urlLonArray->dimension_size(urlLonArray->dim_begin())
94 != urlLonArray->dimension_size(
95 urlLonArray->dim_begin() + 1)) {
97 "The size of the latitude and longitude arrays are not the same");
101 lat.resize(urlLatArray->length());
102 urlLatArray->value(&lat[0]);
103 lon.resize(urlLonArray->length());
104 urlLonArray->value(&lon[0]);
106 VERBOSE(cerr <<
"\tsize of lat array: " << lat.size() << endl);
107 VERBOSE(cerr <<
"\tsize of lon array: " << lon.size() << endl);
109 coords indexVals = coords();
113 vector<float>::iterator j_begin = lon.begin();
115 for (vector<float>::iterator i = lat.begin(), e = lat.end(), j =
116 lon.begin(); i != e; ++i, ++j) {
118 int offset = j - j_begin;
119 indexVals.x = offset / (size_x - 1);
120 indexVals.y = offset % (size_x - 1);
124 indexArray.push_back(indexVals);
127 }
catch (libdap::Error &e) {
128 cerr <<
"ERROR: " << e.get_error_message() << endl;
140 vector<uint64> calculateStareIndex(
const float64 level,
const float64 buildlevel, vector<coords> latlonVals) {
144 STARE index(level,buildlevel);
146 vector<uint64> stareVals;
149 std::unordered_map<float, struct coords> indexMap;
151 for (vector<coords>::iterator i = latlonVals.begin(); i != latlonVals.end(); ++i) {
152 stareIndex = index.ValueFromLatLonDegrees(i->lat, i->lon, level);
153 stareVals.push_back(stareIndex);
156 indexMap[stareIndex] = *i;
165 void writeHDF5(
const string &filename, vector<coords> coordVals, vector<uint64> stareVals) {
166 hid_t file, datasetX, datasetY, datasetLat, datasetLon, datasetStare;
172 vector<float> latArray;
173 vector<float> lonArray;
175 for (vector<coords>::iterator i = coordVals.begin(); i != coordVals.end(); ++i) {
176 xArray.push_back(i->x);
177 yArray.push_back(i->y);
178 latArray.push_back(i->lat);
179 lonArray.push_back(i->lon);
184 hsize_t arrayLength[1] = { xArray.size() };
189 hid_t fapl = H5Pcreate (H5P_FILE_ACCESS);
190 #ifdef H5F_LIBVER_V18
191 H5Pset_libver_bounds (fapl, H5F_LIBVER_EARLIEST, H5F_LIBVER_V18);
193 H5Pset_libver_bounds (fapl, H5F_LIBVER_EARLIEST, H5F_LIBVER_LATEST);
197 file = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, fapl);
201 dataspace = H5Screate_simple(1 , arrayLength, NULL);
213 VERBOSE(cerr <<
"\nCreating datatypes: x, y, lat, lon, stareIndex -> ");
215 dataTypes = H5Tcreate (H5T_COMPOUND,
sizeof(coords));
216 H5Tinsert(dataTypes,
"x", HOFFSET(coords, x), H5T_NATIVE_INT);
217 H5Tinsert(dataTypes,
"y", HOFFSET(coords, y), H5T_NATIVE_INT);
218 H5Tinsert(dataTypes,
"lat", HOFFSET(coords, lat), H5T_NATIVE_FLOAT);
219 H5Tinsert(dataTypes,
"lon", HOFFSET(coords, lon), H5T_NATIVE_FLOAT);
220 H5Tinsert(dataTypes,
"stareIndex", HOFFSET(coords, stareIndex), H5T_NATIVE_INT);
225 const char *datasetName =
"StareData";
226 VERBOSE(cerr <<
"Creating dataset: " << datasetName <<
" -> ");
227 dataset = H5Dcreate2(file, datasetName, dataTypes, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
229 VERBOSE(cerr <<
"Writing data to dataset" << endl);
230 H5Dwrite(dataset, dataTypes, H5S_ALL, H5S_ALL, H5P_DEFAULT, &keyVals[0]);
236 const char *datasetName =
"X";
237 VERBOSE(cerr <<
"Creating dataset: " << datasetName <<
" -> ");
238 datasetX = H5Dcreate2(file, datasetName, H5T_NATIVE_INT, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
241 VERBOSE(cerr <<
"Creating dataset: " << datasetName <<
" -> ");
242 datasetY = H5Dcreate2(file, datasetName, H5T_NATIVE_INT, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
244 datasetName =
"Latitude";
245 VERBOSE(cerr <<
"Creating dataset: " << datasetName <<
" -> ");
246 datasetLat = H5Dcreate2(file, datasetName, H5T_NATIVE_FLOAT, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
248 datasetName =
"Longitude";
249 VERBOSE(cerr <<
"Creating dataset: " << datasetName <<
" -> ");
250 datasetLon = H5Dcreate2(file, datasetName, H5T_NATIVE_FLOAT, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
252 datasetName =
"Stare Index";
253 VERBOSE(cerr <<
"Creating dataset: " << datasetName <<
" -> ");
254 datasetStare = H5Dcreate2(file, datasetName, H5T_NATIVE_INT64, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
259 VERBOSE(cerr <<
"Writing data to dataset" << endl);
260 H5Dwrite(datasetX, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &xArray[0]);
261 H5Dwrite(datasetY, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &yArray[0]);
262 H5Dwrite(datasetLat, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &latArray[0]);
263 H5Dwrite(datasetLon, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &lonArray[0]);
264 H5Dwrite(datasetStare, H5T_NATIVE_INT64, H5S_ALL, H5S_ALL, H5P_DEFAULT, &stareVals[0]);
272 VERBOSE(cerr <<
"Data written to file: " << filename << endl);
282 int main(
int argc,
char *argv[]) {
286 while ((c = getopt(argc, argv,
"hvo:")) != -1) {
289 cerr <<
"\nbuild_sidecar [options] <filename> latitude-name longitude-name level buildlevel\n\n";
290 cerr <<
"-o output file: \tOutput the STARE data to the given output file\n\n";
291 cerr <<
"-v verbose\n" << endl;
304 string dataUrl = argv[argc-5];
305 string latName = argv[argc-4];
306 string lonName = argv[argc-3];
307 float level = atof(argv[argc-2]);
308 float build_level = atof(argv[argc-1]);
317 vector<coords> coordResults = readUrl(dataUrl, latName, lonName);
319 vector<uint64> stareResults = calculateStareIndex(level, build_level, coordResults);
325 vector<coords> keyVals;
326 findLatLon(dataUrl, level, build_level, keyVals);
329 if (newName.empty()) {
333 size_t granulePos = dataUrl.find_last_of(
"/");
334 string granuleName = dataUrl.substr(granulePos + 1);
335 size_t findDot = granuleName.find_last_of(
".");
336 newName = granuleName.substr(0, findDot) +
"_sidecar.h5";
339 writeHDF5(newName, coordResults, stareResults);
341 catch(libdap::Error &e) {
342 cerr <<
"Error: " << e.get_error_message() << endl;