bes  Updated for version 3.20.5
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 "SpatialIndex.h"
14 #include "SpatialInterface.h"
15 
16 #include <D4Connect.h>
17 #include <Connect.h>
18 #include <Response.h>
19 #include <Array.h>
20 
21 #include <Error.h>
22 
23 #include <unordered_map>
24 #include <array>
25 #include <memory>
26 #include "hdf5.h"
27 
28 using namespace std;
29 
30 static bool verbose = false;
31 #define VERBOSE(x) do { if (verbose) x; } while(false)
32 
33 //Struct to store the iterator and lat/lon values that match a specific STARE index
34 struct coords {
35  int x;
36  int y;
37 
38  float64 lat;
39  float64 lon;
40 
41  uint64 stareIndex;
42 };
43 
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) {
53  //Create an htmInterface that will be used to get the STARE index
54  htmInterface htm(level, buildlevel);
55  const SpatialIndex &index = htm.index();
56 
57  auto_ptr < libdap::Connect > url(new libdap::Connect(dataUrl));
58 
59  std::vector<float> lat;
60  std::vector<float> lon;
61 
62  try {
63  libdap::BaseTypeFactory factory;
64  libdap::DataDDS dds(&factory);
65 
66  VERBOSE(cerr << "\n\n\tRequesting data from " << dataUrl << endl);
67 
68  //Makes sure only the Latitude and Longitude variables are requested
69  url->request_data(dds, "Latitude,Longitude");
70 
71  //Create separate libdap arrays to store the lat and lon arrays individually
72  libdap::Array *urlLatArray = dynamic_cast<libdap::Array *>(dds.var(
73  "Latitude"));
74  libdap::Array *urlLonArray = dynamic_cast<libdap::Array *>(dds.var(
75  "Longitude"));
76 
77  // ----Error checking---- //
78  if (urlLatArray == 0 || urlLonArray == 0) {
79  throw libdap::Error("Expected both lat and lon arrays");
80  }
81 
82  unsigned int dims = urlLatArray->dimensions();
83 
84  if (dims != 2) {
85  throw libdap::Error("Incorrect latitude dimensions");
86  }
87 
88  if (dims != urlLonArray->dimensions()) {
89  throw libdap::Error("Incorrect longitude dimensions");
90  }
91 
92  int size_y = urlLatArray->dimension_size(urlLatArray->dim_begin());
93  int size_x = urlLatArray->dimension_size(urlLatArray->dim_begin() + 1);
94 
95  if (size_y != urlLonArray->dimension_size(urlLonArray->dim_begin())
96  || size_x
97  != urlLonArray->dimension_size(
98  urlLonArray->dim_begin() + 1)) {
99  throw libdap::Error(
100  "The size of the latitude and longitude arrays are not the same");
101  }
102 
103  //Initialize the arrays with the correct length for lat and lon
104  lat.resize(urlLatArray->length());
105  urlLatArray->value(&lat[0]);
106  lon.resize(urlLonArray->length());
107  urlLonArray->value(&lon[0]);
108 
109  VERBOSE(cerr << "\tsize of lat array: " << lat.size() << endl);
110  VERBOSE(cerr << "\tsize of lon array: " << lon.size() << endl);
111 
112  coords indexVals = coords();
113  std::unordered_map<float, struct coords> indexMap;
114 
115 #if 0
116  //Taken out because CF doesn't support compound types,
117  // so each variable will need to be stored in its own array
118  // -kln 5/17/19
119 
120  //Array to store the key and values of the indexMap that will be used to write the hdf5 file
121  keyVals.resize(size_x * size_y);
122 #endif
123 
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());
129 
130  //Declare the beginning of the vectors here rather than inside the loop to save compute time
131  //vector<float>::iterator i_begin = lat.begin(); FIXME: Probably not needed since we just use a single dimension array
132  vector<float>::iterator j_begin = lon.begin();
133 
134  int arrayLoc = 0;
135 
136  VERBOSE (cerr << "\nCalculating the STARE indices" << endl);
137 
138  //Make a separate iterator to point at the lat vector and the lon vector inside the same loop.
139  // Then, loop through the lat and lon arrays at the same time
140  for (vector<float>::iterator i = lat.begin(), e = lat.end(), j =
141  lon.begin(); i != e; ++i, ++j) {
142  //Use an offset since we are using a 1D array and treating it like a 2D array
143  int offset = j - j_begin;
144  indexVals.x = offset / (size_x - 1);//Get the current index for the lon
145  indexVals.y = offset % (size_x - 1);//Get the current index for the lat
146  indexVals.lat = *i;
147  indexVals.lon = *j;
148  indexVals.stareIndex = index.idByLatLon(*i, *j);//Use the lat/lon values to calculate the Stare index
149 
150  //Map the STARE index value to the x,y indices
151  indexMap[indexVals.stareIndex] = indexVals;
152 
153  //Store the values calculated for indexVals and store them in their arrays to be used in the hdf5 file
154  xArray[arrayLoc] = indexVals.x;
155  yArray[arrayLoc] = indexVals.y;
156  latArray[arrayLoc] = *i;
157  lonArray[arrayLoc] = *j;
158  stareArray[arrayLoc] = indexVals.stareIndex;
159 
160 #if 0
161  //Taken out because CF doesn't support compound types,
162  // so each variable will need to be stored in its own array
163  // -kln 5/17/19
164 
165  //Store the same previous values inside the coords vector for use in the hdf5 file
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;
171 #endif
172 
173  arrayLoc++;
174  }
175 
176  } catch (libdap::Error &e) {
177  cerr << "ERROR: " << e.get_error_message() << endl;
178  }
179 }
180 
181 
182 /*************
183  * HDF5 Stuff *
184  *************/
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;
187  hid_t dataspace; /* handle */
188 
189  //Used to store the the size of the array.
190  // In other cases where the array is more than 1 dimension each dimension's length would be stored in this array.
191  hsize_t arrayLength[1] = { xArray.size() };
192 
193  //Allows us to use hdf5 1.10 generated files with hdf5 1.8
194  // !---Must be removed if a feature from 1.10 is required to use---!
195  // -kln 5/16/19
196  hid_t fapl = H5Pcreate (H5P_FILE_ACCESS);
197 #ifdef H5F_LIBVER_V18
198  H5Pset_libver_bounds (fapl, H5F_LIBVER_EARLIEST, H5F_LIBVER_V18);
199 #else
200  H5Pset_libver_bounds (fapl, H5F_LIBVER_EARLIEST, H5F_LIBVER_LATEST);
201 #endif
202 
203  //H5Fcreate returns a file id that will be saved in variable "file"
204  file = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, fapl);
205 
206  //The rank is used to determine the dimensions for the dataspace
207  // Since we only use a one dimension array we can always assume the Rank should be 1
208  dataspace = H5Screate_simple(1 /*RANK*/, arrayLength, NULL);
209 
210 #if 0
211  //Taken out because CF doesn't support compound types,
212  // so each variable will need to be stored in its own array
213  // -kln 5/17/19
214 
215  /*
216  * Create the memory datatype.
217  * Because the "coords" struct has ints and floats we need to make sure we put
218  * the correct offset onto memory for each data type
219  */
220  VERBOSE(cerr << "\nCreating datatypes: x, y, lat, lon, stareIndex -> ");
221 
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);
228 
229  /*
230  * Create the dataset
231  */
232  const char *datasetName = "StareData";
233  VERBOSE(cerr << "Creating dataset: " << datasetName << " -> ");
234  dataset = H5Dcreate2(file, datasetName, dataTypes, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
235 
236  VERBOSE(cerr << "Writing data to dataset" << endl);
237  H5Dwrite(dataset, dataTypes, H5S_ALL, H5S_ALL, H5P_DEFAULT, &keyVals[0]);
238 #endif
239 
240  /*
241  * Create the datasets
242  */
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);
246 
247  datasetName = "Y";
248  VERBOSE(cerr << "Creating dataset: " << datasetName << " -> ");
249  datasetY = H5Dcreate2(file, datasetName, H5T_NATIVE_INT, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
250 
251  datasetName = "Latitude";
252  VERBOSE(cerr << "Creating dataset: " << datasetName << " -> ");
253  datasetLat = H5Dcreate2(file, datasetName, H5T_NATIVE_FLOAT, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
254 
255  datasetName = "Longitude";
256  VERBOSE(cerr << "Creating dataset: " << datasetName << " -> ");
257  datasetLon = H5Dcreate2(file, datasetName, H5T_NATIVE_FLOAT, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
258 
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);
262 
263  /*
264  * Write the data to the datasets
265  */
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]);
272 
273  /*
274  * Close/release resources.
275  */
276  H5Sclose(dataspace);
277  H5Fclose(file);
278 
279  VERBOSE(cerr << "Data written to file: " << filename << endl);
280 }
281 
282 
289 int main(int argc, char *argv[]) {
290  int c;
291  string newName;
292 
293  while ((c = getopt(argc, argv, "hvo:")) != -1) {
294  switch (c) {
295  case 'h':
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;
299  exit(1);
300  break;
301  case 'o':
302  newName = optarg;
303  break;
304  case 'v':
305  verbose = true;
306  break;
307  }
308  }
309 
310  //Argument values
311  string dataUrl = argv[argc-3];
312  float level = atof(argv[argc-2]);
313  float build_level = atof(argv[argc-1]);
314 
315  /*if (argc != 4) {
316  cerr << "Expected 3 arguments" << endl;
317  exit(EXIT_FAILURE);
318  }*/
319 
320  try {
321  //Need to store each value in its own array
322  vector<int> xVals;
323  vector<int> yVals;
324  vector<float> latVals;
325  vector<float> lonVals;
326  vector<uint64> stareVals;
327 
328  findLatLon(dataUrl, level, build_level, xVals, yVals, latVals, lonVals, stareVals);
329 #if 0
330  //Taken out because CF doesn't support compound types,
331  // so each variable will need to be stored in its own array
332  // -kln 5/17/19
333 
334  vector<coords> keyVals;
335  findLatLon(dataUrl, level, build_level, keyVals);
336 #endif
337 
338  if (newName.empty()) {
339  //Locate the granule name inside the provided url.
340  // Once the granule is found, add ".h5" to the granule name
341  // Rename the new H5 file to be <granulename.h5>
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";
346  }
347 
348  writeHDF5(newName, xVals, yVals, latVals, lonVals, stareVals);
349  }
350  catch(libdap::Error &e) {
351  cerr << "Error: " << e.get_error_message() << endl;
352  exit(EXIT_FAILURE);
353  }
354 
355  exit(EXIT_SUCCESS);
356 
357 }