bes  Updated for version 3.20.6
stareMapping.cc
1 /*********************************************************************
2  * stareMapping.cc *
3  * *
4  * Created on: Mar 7, 2019 *
5  * *
6  * Purpose: Index the STARE value that corresponds to a given *
7  * iterator from an array, as well as its lat/lon *
8  * value. *
9  * Author: Kodi Neumiller, kneumiller@opendap.org *
10  * *
11  ********************************************************************/
12 
13 #include "STARE.h"
14 
15 #include <D4Connect.h>
16 #include <Connect.h>
17 #include <Response.h>
18 #include <Array.h>
19 
20 #include <Error.h>
21 
22 #include <unordered_map>
23 #include <array>
24 #include <memory>
25 #include "hdf5.h"
26 
27 using namespace std;
28 
29 static bool verbose = false;
30 #define VERBOSE(x) do { if (verbose) x; } while(false)
31 
32 //Struct to store the iterator and lat/lon values that match a specific STARE index
33 struct coords {
34  int x;
35  int y;
36 
37  float64 lat;
38  float64 lon;
39 
40  uint64 stareIndex;
41 };
42 
49 vector<coords> readUrl(string dataUrl, string latName, string lonName) {
50  auto_ptr < libdap::Connect > url(new libdap::Connect(dataUrl));
51 
52  string latlonName = latName + "," + lonName;
53 
54  std::vector<float> lat;
55  std::vector<float> lon;
56 
57  vector<coords> indexArray;
58 
59  try {
60  libdap::BaseTypeFactory factory;
61  libdap::DataDDS dds(&factory);
62 
63  VERBOSE(cerr << "\n\n\tRequesting data from " << dataUrl << endl);
64 
65  //Makes sure only the Latitude and Longitude variables are requested
66  url->request_data(dds, latlonName);
67 
68  //Create separate libdap arrays to store the lat and lon arrays individually
69  libdap::Array *urlLatArray = dynamic_cast<libdap::Array *>(dds.var(
70  latName));
71  libdap::Array *urlLonArray = dynamic_cast<libdap::Array *>(dds.var(
72  lonName));
73 
74  // ----Error checking---- //
75  if (urlLatArray == 0 || urlLonArray == 0) {
76  throw libdap::Error("Expected both lat and lon arrays");
77  }
78 
79  unsigned int dims = urlLatArray->dimensions();
80 
81  if (dims != 2) {
82  throw libdap::Error("Incorrect latitude dimensions");
83  }
84 
85  if (dims != urlLonArray->dimensions()) {
86  throw libdap::Error("Incorrect longitude dimensions");
87  }
88 
89  int size_y = urlLatArray->dimension_size(urlLatArray->dim_begin());
90  int size_x = urlLatArray->dimension_size(urlLatArray->dim_begin() + 1);
91 
92  if (size_y != urlLonArray->dimension_size(urlLonArray->dim_begin())
93  || size_x
94  != urlLonArray->dimension_size(
95  urlLonArray->dim_begin() + 1)) {
96  throw libdap::Error(
97  "The size of the latitude and longitude arrays are not the same");
98  }
99 
100  //Initialize the arrays with the correct length for lat and lon
101  lat.resize(urlLatArray->length());
102  urlLatArray->value(&lat[0]);
103  lon.resize(urlLonArray->length());
104  urlLonArray->value(&lon[0]);
105 
106  VERBOSE(cerr << "\tsize of lat array: " << lat.size() << endl);
107  VERBOSE(cerr << "\tsize of lon array: " << lon.size() << endl);
108 
109  coords indexVals = coords();
110 
111  //Declare the beginning of the vectors here rather than inside the loop to save compute time
112  //vector<float>::iterator i_begin = lat.begin(); FIXME: Probably not needed since we just use a single dimension array
113  vector<float>::iterator j_begin = lon.begin();
114 
115  for (vector<float>::iterator i = lat.begin(), e = lat.end(), j =
116  lon.begin(); i != e; ++i, ++j) {
117  //Use an offset since we are using a 1D array and treating it like a 2D array
118  int offset = j - j_begin;
119  indexVals.x = offset / (size_x - 1);//Get the current index for the lon
120  indexVals.y = offset % (size_x - 1);//Get the current index for the lat
121  indexVals.lat = *i;
122  indexVals.lon = *j;
123 
124  indexArray.push_back(indexVals);
125  }
126 
127  } catch (libdap::Error &e) {
128  cerr << "ERROR: " << e.get_error_message() << endl;
129  }
130 
131  return indexArray;
132 }
133 
140 vector<uint64> calculateStareIndex(const float64 level, const float64 buildlevel, vector<coords> latlonVals) {
141  //Create an htmInterface that will be used to get the STARE index
142  //htmInterface htm(level, buildlevel);
143  //const SpatialIndex &index = htm.index();
144  STARE index(level,buildlevel);
145 
146  vector<uint64> stareVals;
147  uint64 stareIndex;
148 
149  std::unordered_map<float, struct coords> indexMap;
150 
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);
154 
155  //Map the stare values to its corresponding lat/lon and the iterators x/y
156  indexMap[stareIndex] = *i;
157  }
158 
159  return stareVals;
160 }
161 
162 /*************
163  * HDF5 Stuff *
164  *************/
165 void writeHDF5(const string &filename, vector<coords> coordVals, vector<uint64> stareVals) {
166  hid_t file, datasetX, datasetY, datasetLat, datasetLon, datasetStare;
167  hid_t dataspace; /* handle */
168 
169  //Need to store each value in its own array
170  vector<int> xArray;
171  vector<int> yArray;
172  vector<float> latArray;
173  vector<float> lonArray;
174 
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);
180  }
181 
182  //Used to store the the size of the array.
183  // In other cases where the array is more than 1 dimension each dimension's length would be stored in this array.
184  hsize_t arrayLength[1] = { xArray.size() };
185 
186  //Allows us to use hdf5 1.10 generated files with hdf5 1.8
187  // !---Must be removed if a feature from 1.10 is required to use---!
188  // -kln 5/16/19
189  hid_t fapl = H5Pcreate (H5P_FILE_ACCESS);
190 #ifdef H5F_LIBVER_V18
191  H5Pset_libver_bounds (fapl, H5F_LIBVER_EARLIEST, H5F_LIBVER_V18);
192 #else
193  H5Pset_libver_bounds (fapl, H5F_LIBVER_EARLIEST, H5F_LIBVER_LATEST);
194 #endif
195 
196  //H5Fcreate returns a file id that will be saved in variable "file"
197  file = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, fapl);
198 
199  //The rank is used to determine the dimensions for the dataspace
200  // Since we only use a one dimension array we can always assume the Rank should be 1
201  dataspace = H5Screate_simple(1 /*RANK*/, arrayLength, NULL);
202 
203 #if 0
204  //Taken out because CF doesn't support compound types,
205  // so each variable will need to be stored in its own array
206  // -kln 5/17/19
207 
208  /*
209  * Create the memory datatype.
210  * Because the "coords" struct has ints and floats we need to make sure we put
211  * the correct offset onto memory for each data type
212  */
213  VERBOSE(cerr << "\nCreating datatypes: x, y, lat, lon, stareIndex -> ");
214 
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);
221 
222  /*
223  * Create the dataset
224  */
225  const char *datasetName = "StareData";
226  VERBOSE(cerr << "Creating dataset: " << datasetName << " -> ");
227  dataset = H5Dcreate2(file, datasetName, dataTypes, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
228 
229  VERBOSE(cerr << "Writing data to dataset" << endl);
230  H5Dwrite(dataset, dataTypes, H5S_ALL, H5S_ALL, H5P_DEFAULT, &keyVals[0]);
231 #endif
232 
233  /*
234  * Create the datasets
235  */
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);
239 
240  datasetName = "Y";
241  VERBOSE(cerr << "Creating dataset: " << datasetName << " -> ");
242  datasetY = H5Dcreate2(file, datasetName, H5T_NATIVE_INT, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
243 
244  datasetName = "Latitude";
245  VERBOSE(cerr << "Creating dataset: " << datasetName << " -> ");
246  datasetLat = H5Dcreate2(file, datasetName, H5T_NATIVE_FLOAT, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
247 
248  datasetName = "Longitude";
249  VERBOSE(cerr << "Creating dataset: " << datasetName << " -> ");
250  datasetLon = H5Dcreate2(file, datasetName, H5T_NATIVE_FLOAT, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
251 
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);
255 
256  /*
257  * Write the data to the datasets
258  */
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]);
265 
266  /*
267  * Close/release resources.
268  */
269  H5Sclose(dataspace);
270  H5Fclose(file);
271 
272  VERBOSE(cerr << "Data written to file: " << filename << endl);
273 }
274 
275 
282 int main(int argc, char *argv[]) {
283  int c;
284  string newName;
285 
286  while ((c = getopt(argc, argv, "hvo:")) != -1) {
287  switch (c) {
288  case 'h':
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;
292  exit(1);
293  break;
294  case 'o':
295  newName = optarg;
296  break;
297  case 'v':
298  verbose = true;
299  break;
300  }
301  }
302 
303  //Argument values
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]);
309 
310  /*if (argc != 4) {
311  cerr << "Expected 3 arguments" << endl;
312  exit(EXIT_FAILURE);
313  }*/
314 
315  try {
316 
317  vector<coords> coordResults = readUrl(dataUrl, latName, lonName);
318 
319  vector<uint64> stareResults = calculateStareIndex(level, build_level, coordResults);
320 #if 0
321  //Taken out because CF doesn't support compound types,
322  // so each variable will need to be stored in its own array
323  // -kln 5/17/19
324 
325  vector<coords> keyVals;
326  findLatLon(dataUrl, level, build_level, keyVals);
327 #endif
328 
329  if (newName.empty()) {
330  //Locate the granule name inside the provided url.
331  // Once the granule is found, add ".h5" to the granule name
332  // Rename the new H5 file to be <granulename.h5>
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";
337  }
338 
339  writeHDF5(newName, coordResults, stareResults);
340  }
341  catch(libdap::Error &e) {
342  cerr << "Error: " << e.get_error_message() << endl;
343  exit(EXIT_FAILURE);
344  }
345 
346  exit(EXIT_SUCCESS);
347 
348 }