bes  Updated for version 3.20.6
HDFEOS2ArrayGridGeoField.h
1 // This file is part of the hdf4 data handler for the OPeNDAP data server.
3 // It retrieves the latitude and longitude of the HDF-EOS2 Grid
4 // There are two typical cases:
5 // read the lat/lon from the file and calculate lat/lon using the EOS2 library.
6 // Several variations are also handled:
7 // 1. For geographic and Cylinderic Equal Area projections, condense 2-D to 1-D.
8 // 2. For some files, the longitude is within 0-360 range instead of -180 - 180 range.
9 // We need to convert 0-360 to -180-180.
10 // 3. Some files have fillvalues in the lat. and lon. for the geographic projection.
11 // 4. Several MODIS files don't have the correct parameters inside StructMetadata.
12 // We can obtain the starting point, the step and replace the fill value.
13 // Authors: MuQun Yang <myang6@hdfgroup.org> Choonghwan Lee
14 // Copyright (c) 2009-2012 The HDF Group
16 #ifdef USE_HDFEOS2_LIB
17 #ifndef HDFEOS2ARRAY_GRIDGEOFIELD_H
18 #define HDFEOS2ARRAY_GRIDGEOFIELD_H
19 
20 #include <cmath>
21 
22 #include <Array.h>
23 
24 #include "mfhdf.h"
25 #include "hdf.h"
26 #include "HdfEosDef.h"
27 
28 
29 class HDFEOS2ArrayGridGeoField:public libdap::Array
30 {
31  public:
32  HDFEOS2ArrayGridGeoField (int rank, int fieldtype, bool llflag, bool ydimmajor, bool condenseddim, bool speciallon, int specialformat, /*short field_cache,*/const std::string &filename, const int gridfd, const std::string & gridname, const std::string & fieldname,const string & n = "", libdap::BaseType * v = 0):
33  libdap::Array (n, v),
34  rank (rank),
35  fieldtype (fieldtype),
36  llflag (llflag),
37  ydimmajor (ydimmajor),
38  condenseddim (condenseddim),
39  speciallon (speciallon),
40  specialformat (specialformat),
41  /*field_cache(field_cache),*/
42  filename(filename),
43  gridfd(gridfd),
44  gridname (gridname),
45  fieldname (fieldname)
46  {
47  }
48  virtual ~ HDFEOS2ArrayGridGeoField ()
49  {
50  }
51  int format_constraint (int *cor, int *step, int *edg);
52 
53  libdap::BaseType *ptr_duplicate ()
54  {
55  return new HDFEOS2ArrayGridGeoField (*this);
56  }
57 
58  virtual bool read ();
59 
60  private:
61 
62  // Field array rank
63  int rank;
64 
65  // Distinguish coordinate variables from general variables.
66  // For fieldtype values:
67  // 0 the field is a general field
68  // 1 the field is latitude.
69  // 2 the field is longtitude.
70  // 3 the field is a coordinate variable defined as level.
71  // 4 the field is an inserted natural number.
72  // 5 the field is time.
73  int fieldtype;
74 
75  // The flag to indicate if lat/lon is an existing field in the file or needs to be calculated.
76  bool llflag;
77 
78  // Flag to check if this lat/lon field is YDim major(YDim,XDim). This is necessary to use GDij2ll
79  bool ydimmajor;
80 
81  // Flag to check if this 2-D lat/lon can be condensed to 1-D lat/lon
82  bool condenseddim;
83 
84  // Flag to check if this file's longitude needs to be handled specially.
85  // Note: longitude values range from 0 to 360 for some fields. We need to map the values to -180 to 180.
86  bool speciallon;
87 
88  // Latitude and longitude values of some HDF-EOS2 grids need to be handled in special ways.
89  // There are four cases that we need to calculate lat/lon differently.
90  // This number is used to distinguish them.
91  // 1) specialformat = 1
92  // Projection: Geographic
93  // upleft and lowright coordinates don't follow EOS's DDDMMMSSS conventions.
94  // Instead, they simply represent lat/lon values as -180.0 or -90.0.
95  // Products: mostly MODIS MCD Grid
96 
97  // 2) specialformat = 2
98  // Projection: Geographic
99  // upleft and lowright coordinates don't follow EOS's DDDMMMSSS conventions.
100  // Instead, their upleft or lowright are simply represented as default(-1).
101  // Products: mostly TRMM CERES Grid
102 
103  // 3) specialformat = 3
104  // One MOD13C2 doesn't provide project code
105  // The upperleft and lowerright coordinates are all -1
106  // We have to calculate lat/lon by ourselves.
107  // Since it doesn't provide the project code, we double check their information
108  // and find that it covers the whole globe with 0.05 degree resolution.
109  // Lat. is from 90 to -90 and Lon is from -180 to 180.
110 
111  // 4) specialformat = 4
112  // Projection: Space Oblique Mercator(SOM)
113  // The lat/lon needs to be handled differently for the SOM projection
114  // Products: MISR
115  int specialformat;
116 
117  //short field_cache;
118 
119  // Temp here: HDF-EOS2 file name
120  std::string filename;
121 
122  int gridfd;
123 
124  // HDF-EOS2 grid name
125  std::string gridname;
126 
127  // HDF-EOS2 field name
128  std::string fieldname;
129  // Calculate Lat and Lon based on HDF-EOS2 library.
130  void CalculateLatLon (int32 gridid, int fieldtype, int specialformat, float64 * outlatlon, float64* latlon_all, int32 * offset, int32 * count, int32 * step, int nelms,bool write_latlon_cache);
131 
132  // Calculate Special Latitude and Longitude.
133  //One MOD13C2 file doesn't provide projection code
134  // The upperleft and lowerright coordinates are all -1
135  // We have to calculate lat/lon by ourselves.
136  // Since it doesn't provide the project code, we double check their information
137  // and find that it covers the whole globe with 0.05 degree resolution.
138  // Lat. is from 90 to -90 and Lon is from -180 to 180.
139  void CalculateSpeLatLon (int32 gridid, int fieldtype, float64 * outlatlon, int32 * offset, int32 * count, int32 * step);
140 
141  // Calculate Latitude and Longtiude for the Geo-projection for very large number of elements per dimension.
142  void CalculateLargeGeoLatLon(int32 gridid, int fieldtype, float64* latlon, float64* latlon_all, int *start, int *count, int *step, int nelms,bool write_latlon_cache);
143  // test for undefined values returned by longitude-latitude calculation
144  bool isundef_lat(double value)
145  {
146  if (std::isinf(value)) return(true);
147  if (std::isnan(value)) return(true);
148  // GCTP_AMAZ returns "1e+51" for values at the opposite poles
149  if(value < -90.0 || value > 90.0)
150  return(true);
151  // This is ok.
152  return(false);
153  } // end bool isundef_lat(double value)
154 
155  bool isundef_lon(double value)
156  {
157  if (std::isinf(value)) return(true);
158  if (std::isnan(value)) return(true);
159  // GCTP_LAMAZ returns "1e+51" for values at the opposite poles
160  if(value < -180.0 || value > 180.0) return(true);
161  return(false);
162  } // end bool isundef_lat(double value)
163 
164  // Given rol, col address in double array of dimension YDim x XDim
165  // return value of nearest neighbor to (row,col) which is not undefined
166  double nearestNeighborLatVal(double *array, int row, int col, int YDim, int XDim)
167  {
168  // test valid row, col address range
169  if(row < 0 || row > YDim || col < 0 || col > XDim)
170  {
171  cerr << "nearestNeighborLatVal("<<row<<", "<<col<<", "<<YDim<<", "<<XDim;
172  cerr <<"): index out of range"<<endl;
173  return(0.0);
174  }
175  // address (0,0)
176  if(row < YDim/2 && col < XDim/2)
177  { /* search by incrementing both row and col */
178  if(!isundef_lat(array[(row+1)*XDim+col])) return(array[(row+1)*XDim+col]);
179  if(!isundef_lat(array[row*XDim+col+1])) return(array[row*XDim+col+1]);
180  if(!isundef_lat(array[(row+1)*XDim+col+1])) return(array[(row+1)*XDim+col+1]);
181  /* recurse on the diagonal */
182  return(nearestNeighborLatVal(array, row+1, col+1, YDim, XDim));
183  }
184  if(row < YDim/2 && col > XDim/2)
185  { /* search by incrementing row and decrementing col */
186  if(!isundef_lat(array[(row+1)*XDim+col])) return(array[(row+1)*XDim+col]);
187  if(!isundef_lat(array[row*XDim+col-1])) return(array[row*XDim+col-1]);
188  if(!isundef_lat(array[(row+1)*XDim+col-1])) return(array[(row+1)*XDim+col-1]);
189  /* recurse on the diagonal */
190  return(nearestNeighborLatVal(array, row+1, col-1, YDim, XDim));
191  }
192  if(row > YDim/2 && col < XDim/2)
193  { /* search by incrementing col and decrementing row */
194  if(!isundef_lat(array[(row-1)*XDim+col])) return(array[(row-1)*XDim+col]);
195  if(!isundef_lat(array[row*XDim+col+1])) return(array[row*XDim+col+1]);
196  if(!isundef_lat(array[(row-1)*XDim+col+1])) return(array[(row-1)*XDim+col+1]);
197  /* recurse on the diagonal */
198  return(nearestNeighborLatVal(array, row-1, col+1, YDim, XDim));
199  }
200  if(row > YDim/2 && col > XDim/2)
201  { /* search by decrementing both row and col */
202  if(!isundef_lat(array[(row-1)*XDim+col])) return(array[(row-1)*XDim+col]);
203  if(!isundef_lat(array[row*XDim+col-1])) return(array[row*XDim+col-1]);
204  if(!isundef_lat(array[(row-1)*XDim+col-1])) return(array[(row-1)*XDim+col-1]);
205  /* recurse on the diagonal */
206  return(nearestNeighborLatVal(array, row-1, col-1, YDim, XDim));
207  }
208  // dummy return, turn off the compiling warning
209  return 0.0;
210  } // end
211 
212  double nearestNeighborLonVal(double *array, int row, int col, int YDim, int XDim)
213  {
214  // test valid row, col address range
215  if(row < 0 || row > YDim || col < 0 || col > XDim)
216  {
217  cerr << "nearestNeighborLonVal("<<row<<", "<<col<<", "<<YDim<<", "<<XDim;
218  cerr <<"): index out of range"<<endl;
219  return(0.0);
220  }
221  // address (0,0)
222  if(row < YDim/2 && col < XDim/2)
223  { /* search by incrementing both row and col */
224  if(!isundef_lon(array[(row+1)*XDim+col])) return(array[(row+1)*XDim+col]);
225  if(!isundef_lon(array[row*XDim+col+1])) return(array[row*XDim+col+1]);
226  if(!isundef_lon(array[(row+1)*XDim+col+1])) return(array[(row+1)*XDim+col+1]);
227  /* recurse on the diagonal */
228  return(nearestNeighborLonVal(array, row+1, col+1, YDim, XDim));
229  }
230  if(row < YDim/2 && col > XDim/2)
231  { /* search by incrementing row and decrementing col */
232  if(!isundef_lon(array[(row+1)*XDim+col])) return(array[(row+1)*XDim+col]);
233  if(!isundef_lon(array[row*XDim+col-1])) return(array[row*XDim+col-1]);
234  if(!isundef_lon(array[(row+1)*XDim+col-1])) return(array[(row+1)*XDim+col-1]);
235  /* recurse on the diagonal */
236  return(nearestNeighborLonVal(array, row+1, col-1, YDim, XDim));
237  }
238  if(row > YDim/2 && col < XDim/2)
239  { /* search by incrementing col and decrementing row */
240  if(!isundef_lon(array[(row-1)*XDim+col])) return(array[(row-1)*XDim+col]);
241  if(!isundef_lon(array[row*XDim+col+1])) return(array[row*XDim+col+1]);
242  if(!isundef_lon(array[(row-1)*XDim+col+1])) return(array[(row-1)*XDim+col+1]);
243  /* recurse on the diagonal */
244  return(nearestNeighborLonVal(array, row-1, col+1, YDim, XDim));
245  }
246  if(row > YDim/2 && col > XDim/2)
247  { /* search by decrementing both row and col */
248  if(!isundef_lon(array[(row-1)*XDim+col])) return(array[(row-1)*XDim+col]);
249  if(!isundef_lon(array[row*XDim+col-1])) return(array[row*XDim+col-1]);
250  if(!isundef_lon(array[(row-1)*XDim+col-1])) return(array[(row-1)*XDim+col-1]);
251  /* recurse on the diagonal */
252  return(nearestNeighborLonVal(array, row-1, col-1, YDim, XDim));
253  }
254 
255  // dummy return, turn off the compiling warning
256  return 0.0;
257  } // end
258 
259  // Calculate Latitude and Longitude for SOM Projection.
260  // since the latitude and longitude of the SOM projection are 3-D, so we need to handle this projection in a special way.
261  // Based on our current understanding, the third dimension size is always 180.
262  // This is according to the MISR Lat/lon calculation document
263  // at http://eosweb.larc.nasa.gov/PRODOCS/misr/DPS/DPS_v50_RevS.pdf
264  void CalculateSOMLatLon(int32, int*, int*, int*, int,const string &, bool);
265 
266  // Calculate Latitude and Longitude for LAMAZ Projection.
267  void CalculateLAMAZLatLon(int32, int, float64*, float64*,int*, int*, int*, bool);
268 
269  // Subsetting the latitude and longitude.
270  template <class T> void LatLon2DSubset (T* outlatlon, int ydim, int xdim, T* latlon, int32 * offset, int32 * count, int32 * step);
271 
272  // Handle latitude and longitude when having fill value for geographic projection
273  //template <class T> void HandleFillLatLon(T* total_latlon, T* latlon,bool ydimmajor,
274  template <class T> void HandleFillLatLon(vector<T> total_latlon, T* latlon,bool ydimmajor, int fieldtype, int32 xdim , int32 ydim, int32* offset, int32* count, int32* step, int fv);
275 
276  // Corrected Latitude and longitude when the lat/lon has fill value case.
277  template < class T > bool CorLatLon (T * latlon, int fieldtype, int elms, int fv);
278 
279  // Converted longitude from 0-360 to -180-180.
280  template < class T > void CorSpeLon (T * lon, int xdim);
281 
282  // Lat and Lon for GEO and CEA projections need to be condensed from 2-D to 1-D.
283  // This function does this.
284  void getCorrectSubset (int *offset, int *count, int *step, int32 * offset32, int32 * count32, int32 * step32, bool condenseddim, bool ydimmajor, int fieldtype, int rank);
285 
286  // Helper function to handle the case that lat. and lon. contain fill value.
287  template < class T > int findfirstfv (T * array, int start, int end, int fillvalue);
288 
289 };
290 #endif
291 #endif
libdap
Definition: BESDapFunctionResponseCache.h:35