bes  Updated for version 3.20.6
FONgGrid.cc
1 // FONgGrid.cc
2 
3 // This file is part of BES GDAL File Out Module
4 
5 // Copyright (c) 2012 OPeNDAP, Inc.
6 // Author: James Gallagher <jgallagher@opendap.org>
7 //
8 // This library is free software; you can redistribute it and/or
9 // modify it under the terms of the GNU Lesser General Public
10 // License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21 //
22 // You can contact University Corporation for Atmospheric Research at
23 // 3080 Center Green Drive, Boulder, CO 80301
24 
25 #include <algorithm>
26 
27 #include <gdal.h>
28 #include <gdal_priv.h>
29 #include <ogr_spatialref.h>
30 
31 #include <DDS.h>
32 #include <Grid.h>
33 // #include <ce_functions.h>
34 #include <util.h>
35 
36 #include <BESInternalError.h>
37 #include <BESDebug.h>
38 
39 #include "GeoTiffTransmitter.h"
40 #include "FONgTransform.h"
41 #include "FONgGrid.h"
42 
43 using namespace libdap;
44 
49 FONgGrid::FONgGrid(Grid *g): d_grid(g), d_lat(0), d_lon(0)
50 
51 {
52  d_type = dods_grid_c;
53 
54  // Build sets of attribute values for easy searching.
55  // Copied from GeoConstriant in libdap. That class is
56  // abstract and didn't want to modify libdap's ABI for this hack.
57 
58  d_coards_lat_units.insert("degrees_north");
59  d_coards_lat_units.insert("degree_north");
60  d_coards_lat_units.insert("degree_N");
61  d_coards_lat_units.insert("degrees_N");
62 
63  d_coards_lon_units.insert("degrees_east");
64  d_coards_lon_units.insert("degree_east");
65  d_coards_lon_units.insert("degrees_E");
66  d_coards_lon_units.insert("degree_E");
67 
68  d_lat_names.insert("COADSY");
69  d_lat_names.insert("lat");
70  d_lat_names.insert("Lat");
71  d_lat_names.insert("LAT");
72 
73  d_lon_names.insert("COADSX");
74  d_lon_names.insert("lon");
75  d_lon_names.insert("Lon");
76  d_lon_names.insert("LON");
77 }
78 
85 {
86 }
87 
92 class is_prefix {
93 private:
94  string s;
95 public:
96  is_prefix(const string & in) :
97  s(in)
98  {
99  }
100 
101  bool operator()(const string & prefix)
102  {
103  return s.find(prefix) == 0;
104  }
105 };
106 
117 bool FONgGrid::m_lat_unit_or_name_match(const string &var_units, const string &var_name, const string &long_name)
118 {
119  return (long_name == "latitude" || d_coards_lat_units.find(var_units) != d_coards_lat_units.end()
120  || find_if(d_lat_names.begin(), d_lat_names.end(), is_prefix(var_name)) != d_lat_names.end());
121 }
122 
123 bool FONgGrid::m_lon_unit_or_name_match(const string &var_units, const string &var_name, const string &long_name)
124 {
125  return (long_name == "longitude" || d_coards_lon_units.find(var_units) != d_coards_lon_units.end()
126  || find_if(d_lon_names.begin(), d_lon_names.end(), is_prefix(var_name)) != d_lon_names.end());
127 }
128 
145 {
146  Grid::Map_iter m = d_grid->map_begin();
147 
148  // Assume that a Grid is correct and thus has exactly as many maps as its
149  // array part has dimensions. Thus, don't bother to test the Grid's array
150  // dimension iterator for '!= dim_end()'.
151  Array::Dim_iter d = d_grid->get_array()->dim_begin();
152 
153  // The fields d_latitude and d_longitude may be initialized to null or they
154  // may already contain pointers to the maps to use. In the latter case,
155  // skip the heuristics used in this code. However, given that all this
156  // method does is find the lat/lon maps, if they are given in the ctor,
157  // This method will likely not be called at all.
158 
159  while (m != d_grid->map_end() && (!d_lat || !d_lon)) {
160  string units_value = (*m)->get_attr_table().get_attr("units");
161  units_value = remove_quotes(units_value);
162  string long_name = (*m)->get_attr_table().get_attr("long_name");
163  long_name = remove_quotes(long_name);
164  string map_name = (*m)->name();
165 
166  // The 'units' attribute must match exactly; the name only needs to
167  // match a prefix.
168  if (!d_lat && m_lat_unit_or_name_match(units_value, map_name, long_name)) {
169  // Set both d_latitude (a pointer to the real map vector) and
170  // d_lat, a vector of the values represented as doubles. It's easier
171  // to work with d_lat, but it's d_latitude that needs to be set
172  // when constraining the grid. Also, record the grid variable's
173  // dimension iterator so that it's easier to set the Grid's Array
174  // (which also has to be constrained).
175  d_lat = dynamic_cast<Array *>(*m);
176  if (!d_lat) throw InternalErr(__FILE__, __LINE__, "Expected an array.");
177 
178  if (!d_lat->read_p()) d_lat->read();
179  }
180 
181  if (!d_lon && m_lon_unit_or_name_match(units_value, map_name, long_name)) {
182  d_lon = dynamic_cast<Array *>(*m);
183  if (!d_lon) throw InternalErr(__FILE__, __LINE__, "Expected an array.");
184 
185  if (!d_lon->read_p()) d_lon->read();
186  }
187 
188  ++m;
189  ++d;
190  }
191 
192  return d_lat && d_lon;
193 }
194 
204 {
205  BESDEBUG("fong3", "Entering FONgGrid::extract_coordinates" << endl);
206 
207  double *lat = 0, *lon = 0;
208  try {
209  // Find the lat and lon maps for this Grid; throws Error
210  // if the are not found.
212 
213  // The array size
214  t.set_height(d_lat->length());
215  t.set_width(d_lon->length());
216 
217  lat = extract_double_array(d_lat);
218  lon = extract_double_array(d_lon);
219 
220  //max latidude
221  t.set_top(max (lat[0], lat[t.height()-1]));
222  t.set_left(lon[0]);
223  // min latitude
224  t.set_bottom(min (lat[0], lat[t.height()-1]));
225  t.set_right(lon[t.width() - 1]);
226 
227  // Read this from the 'missing_value' or '_FillValue' attributes
228  string missing_value = d_grid->get_attr_table().get_attr("missing_value");
229  if (missing_value.empty()) missing_value = d_grid->get_array()->get_attr_table().get_attr("missing_value");
230  if (missing_value.empty()) missing_value = d_grid->get_attr_table().get_attr("_FillValue");
231  if (missing_value.empty()) missing_value = d_grid->get_array()->get_attr_table().get_attr("_FillValue");
232 
233  BESDEBUG("fong3", "missing_value attribute: " << missing_value << endl);
234 
235  // NB: no_data_type() is 'none' by default
236  if (!missing_value.empty()) {
237  t.set_no_data(missing_value);
238  if (t.no_data() < 0)
239  t.set_no_data_type(FONgTransform::negative);
240  else
241  t.set_no_data_type(FONgTransform::positive);
242  }
243 
244  t.geo_transform_set(true);
245 
246  t.set_num_bands(t.num_bands() + 1);
247  t.push_var(this);
248 
249  delete[] lat;
250  delete[] lon;
251  }
252  catch (Error &e) {
253  delete[] lat;
254  delete[] lon;
255 
256  throw;
257  }
258 }
259 
261 static bool is_spherical(BaseType *btp)
262 {
263  /* crs:grid_mapping_name = "latitude_longitude"
264  crs:semi_major_axis = 6371000.0 ;
265  crs:inverse_flattening = 0 ; */
266 
267  bool gmn = btp->get_attr_table().get_attr("grid_mapping_name") == "latitude_longitude";
268  bool sma = btp->get_attr_table().get_attr("semi_major_axis") == "6371000.0";
269  bool iflat = btp->get_attr_table().get_attr("inverse_flattening") == "0";
270 
271  return gmn && sma && iflat;
272 }
273 
275 static bool is_wgs84(BaseType *btp)
276 {
277  /*crs:grid_mapping_name = "latitude_longitude";
278  crs:longitude_of_prime_meridian = 0.0 ;
279  crs:semi_major_axis = 6378137.0 ;
280  crs:inverse_flattening = 298.257223563 ; */
281 
282  bool gmn = btp->get_attr_table().get_attr("grid_mapping_name") == "latitude_longitude";
283  bool lpm = btp->get_attr_table().get_attr("longitude_of_prime_meridian") == "0.0";
284  bool sma = btp->get_attr_table().get_attr("semi_major_axis") == "6378137.0";
285  bool iflat = btp->get_attr_table().get_attr("inverse_flattening") == "298.257223563";
286 
287  return gmn && lpm && sma && iflat;
288 }
289 
296 string FONgGrid::get_projection(DDS *dds)
297 {
298  // Here's the information about the CF and projections
299  // http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#grid-mappings-and-projections
300  // How this code looks for mapping information: Look for an
301  // attribute named 'grid_mapping' and get it's value. This attribute
302  // with be the name of a variable in the dataset, so get that
303  // variable. Now, look at the attributes of that variable.
304  string mapping_info = d_grid->get_attr_table().get_attr("grid_mapping");
305  if (mapping_info.empty()) mapping_info = d_grid->get_array()->get_attr_table().get_attr("grid_mapping");
306 
307  string WK_GCS = GeoTiffTransmitter::default_gcs;
308 
309  if (!mapping_info.empty()) {
310  // "WGS84": same as "EPSG:4326" but has no dependence on EPSG data files.
311  // "WGS72": same as "EPSG:4322" but has no dependence on EPSG data files.
312  // "NAD27": same as "EPSG:4267" but has no dependence on EPSG data files.
313  // "NAD83": same as "EPSG:4269" but has no dependence on EPSG data files.
314  // "EPSG:n": same as doing an ImportFromEPSG(n).
315 
316  // The mapping info is actually stored as attributes of an Int32 variable.
317  BaseType *btp = dds->var(mapping_info);
318  if (btp && btp->name() == "crs") {
319  if (is_wgs84(btp))
320  WK_GCS = "WGS84";
321  else if (is_spherical(btp)) WK_GCS = "EPSG:4047";
322  }
323  }
324 
325  OGRSpatialReference srs;
326  srs.SetWellKnownGeogCS(WK_GCS.c_str());
327  char *srs_wkt = NULL;
328  srs.exportToWkt(&srs_wkt);
329 
330  string wkt = srs_wkt; // save the result
331 
332  CPLFree(srs_wkt); // free memory alloc'd by GDAL
333 
334  return wkt;
335 }
336 
338 {
339  if (!d_grid->get_array()->read_p()) d_grid->get_array()->read();
340 
341  // This code assumes read() has been called.
342  return extract_double_array(d_grid->get_array());
343 }
344 
FONgGrid::extract_coordinates
virtual void extract_coordinates(FONgTransform &t)
Get the GDAL/OGC WKT projection string.
Definition: FONgGrid.cc:203
FONgGrid::find_lat_lon_maps
bool find_lat_lon_maps()
Definition: FONgGrid.cc:144
libdap
Definition: BESDapFunctionResponseCache.h:35
FONgGrid::get_data
virtual double * get_data()
Get the data values for the band(s). Call must delete.
Definition: FONgGrid.cc:337
FONgGrid::~FONgGrid
virtual ~FONgGrid()
Destructor that cleans up the grid.
Definition: FONgGrid.cc:84
FONgGrid::FONgGrid
FONgGrid(libdap::Grid *g)
Constructor for FONgGrid that takes a DAP Grid.
Definition: FONgGrid.cc:49
FONgTransform
Transformation object that converts an OPeNDAP DataDDS to a GeoTiff file.
Definition: FONgTransform.h:41
Error
FONgGrid::get_projection
string get_projection(libdap::DDS *dds)
Set the projection information For Grids, look for CF information. If it's not present,...
Definition: FONgGrid.cc:296