bes Updated for version 3.20.10
HDFEOS2ArrayGridGeoField.h
1
2// 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 <libdap/Array.h>
23
24#include "mfhdf.h"
25#include "hdf.h"
26#include "HdfEosDef.h"
27
28
29class 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