FONgGrid.cc

// FONgGrid.cc
 
// This file is part of BES GDAL File Out Module
 
// Copyright (c) 2012 OPeNDAP, Inc.
// Author: James Gallagher <jgallagher@opendap.org>
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
//
// You can contact University Corporation for Atmospheric Research at
// 3080 Center Green Drive, Boulder, CO 80301
 
#include <algorithm>
 
#include <gdal.h>
#include <gdal_priv.h>
#include <ogr_spatialref.h>
 
#include <libdap/DDS.h>
#include <libdap/Grid.h>
// #include <ce_functions.h>
#include <libdap/util.h>
 
#include <BESInternalError.h>
#include <BESDebug.h>
 
#include "GeoTiffTransmitter.h"
#include "FONgTransform.h"
#include "FONgGrid.h"
 
using namespace libdap;

FONgGrid::FONgGrid(Grid *g): d_grid(g), d_lat(0), d_lon(0)
 
{
    d_type = dods_grid_c;
 
    // Build sets of attribute values for easy searching.
    // Copied from GeoConstriant in libdap. That class is
    // abstract and didn't want to modify libdap's ABI for this hack.
 
    d_coards_lat_units.insert("degrees_north");
    d_coards_lat_units.insert("degree_north");
    d_coards_lat_units.insert("degree_N");
    d_coards_lat_units.insert("degrees_N");
 
    d_coards_lon_units.insert("degrees_east");
    d_coards_lon_units.insert("degree_east");
    d_coards_lon_units.insert("degrees_E");
    d_coards_lon_units.insert("degree_E");
 
    d_lat_names.insert("COADSY");
    d_lat_names.insert("lat");
    d_lat_names.insert("Lat");
    d_lat_names.insert("LAT");
 
    d_lon_names.insert("COADSX");
    d_lon_names.insert("lon");
    d_lon_names.insert("Lon");
    d_lon_names.insert("LON");
}

FONgGrid::~FONgGrid()
{
}

class is_prefix {
private:
    string s;
public:
    is_prefix(const string & in) :
        s(in)
    {
    }
 
    bool operator()(const string & prefix)
    {
        return s.find(prefix) == 0;
    }
};

bool FONgGrid::m_lat_unit_or_name_match(const string &var_units, const string &var_name, const string &long_name)
{
    return (long_name == "latitude" || d_coards_lat_units.find(var_units) != d_coards_lat_units.end()
        || find_if(d_lat_names.begin(), d_lat_names.end(), is_prefix(var_name)) != d_lat_names.end());
}
 
bool FONgGrid::m_lon_unit_or_name_match(const string &var_units, const string &var_name, const string &long_name)
{
    return (long_name == "longitude" || d_coards_lon_units.find(var_units) != d_coards_lon_units.end()
        || find_if(d_lon_names.begin(), d_lon_names.end(), is_prefix(var_name)) != d_lon_names.end());
}

bool FONgGrid::find_lat_lon_maps()
{
    Grid::Map_iter m = d_grid->map_begin();
 
    // Assume that a Grid is correct and thus has exactly as many maps as its
    // array part has dimensions. Thus, don't bother to test the Grid's array
    // dimension iterator for '!= dim_end()'.
    Array::Dim_iter d = d_grid->get_array()->dim_begin();
 
    // The fields d_latitude and d_longitude may be initialized to null or they
    // may already contain pointers to the maps to use. In the latter case,
    // skip the heuristics used in this code. However, given that all this
    // method does is find the lat/lon maps, if they are given in the ctor,
    // This method will likely not be called at all.
 
    while (m != d_grid->map_end() && (!d_lat || !d_lon)) {
        string units_value = (*m)->get_attr_table().get_attr("units");
        units_value = remove_quotes(units_value);
        string long_name = (*m)->get_attr_table().get_attr("long_name");
        long_name = remove_quotes(long_name);
        string map_name = (*m)->name();
 
        // The 'units' attribute must match exactly; the name only needs to
        // match a prefix.
        if (!d_lat && m_lat_unit_or_name_match(units_value, map_name, long_name)) {
            // Set both d_latitude (a pointer to the real map vector) and
            // d_lat, a vector of the values represented as doubles. It's easier
            // to work with d_lat, but it's d_latitude that needs to be set
            // when constraining the grid. Also, record the grid variable's
            // dimension iterator so that it's easier to set the Grid's Array
            // (which also has to be constrained).
            d_lat = dynamic_cast<Array *>(*m);
            if (!d_lat) throw InternalErr(__FILE__, __LINE__, "Expected an array.");
 
            if (!d_lat->read_p()) d_lat->read();
        }
 
        if (!d_lon && m_lon_unit_or_name_match(units_value, map_name, long_name)) {
            d_lon = dynamic_cast<Array *>(*m);
            if (!d_lon) throw InternalErr(__FILE__, __LINE__, "Expected an array.");
 
            if (!d_lon->read_p()) d_lon->read();
        }
 
        ++m;
        ++d;
    }
 
    return d_lat && d_lon;
}

void FONgGrid::extract_coordinates(FONgTransform &t)
{
    BESDEBUG("fong3", "Entering FONgGrid::extract_coordinates" << endl);
 
    double *lat = 0, *lon = 0;
    try {
        // Find the lat and lon maps for this Grid; throws Error
        // if the are not found.
        find_lat_lon_maps();
 
        // The array size
        t.set_height(d_lat->length());
        t.set_width(d_lon->length());
 
        lat = extract_double_array(d_lat);
        lon = extract_double_array(d_lon);
 
        //max latidude
        t.set_top(max (lat[0], lat[t.height()-1]));
        t.set_left(lon[0]);
        // min latitude
        t.set_bottom(min (lat[0], lat[t.height()-1]));
        t.set_right(lon[t.width() - 1]);
 
        // Read this from the 'missing_value' or '_FillValue' attributes
        string missing_value = d_grid->get_attr_table().get_attr("missing_value");
        if (missing_value.empty()) missing_value = d_grid->get_array()->get_attr_table().get_attr("missing_value");
        if (missing_value.empty()) missing_value = d_grid->get_attr_table().get_attr("_FillValue");
        if (missing_value.empty()) missing_value = d_grid->get_array()->get_attr_table().get_attr("_FillValue");
 
        BESDEBUG("fong3", "missing_value attribute: " << missing_value << endl);
 
        // NB: no_data_type() is 'none' by default
        if (!missing_value.empty()) {
            t.set_no_data(missing_value);
            if (t.no_data() < 0)
                t.set_no_data_type(FONgTransform::negative);
            else
                t.set_no_data_type(FONgTransform::positive);
        }
 
        t.geo_transform_set(true);
 
        t.set_num_bands(t.num_bands() + 1);
        t.push_var(this);
 
        delete[] lat;
        delete[] lon;
    }
    catch (Error &e) {
        delete[] lat;
        delete[] lon;
 
        throw;
    }
}

static bool is_spherical(BaseType *btp)
{
    /* crs:grid_mapping_name = "latitude_longitude"
     crs:semi_major_axis = 6371000.0 ;
     crs:inverse_flattening = 0 ; */
 
    bool gmn = btp->get_attr_table().get_attr("grid_mapping_name") == "latitude_longitude";
    bool sma = btp->get_attr_table().get_attr("semi_major_axis") == "6371000.0";
    bool iflat = btp->get_attr_table().get_attr("inverse_flattening") == "0";
 
    return gmn && sma && iflat;
}

static bool is_wgs84(BaseType *btp)
{
    /*crs:grid_mapping_name = "latitude_longitude";
     crs:longitude_of_prime_meridian = 0.0 ;
     crs:semi_major_axis = 6378137.0 ;
     crs:inverse_flattening = 298.257223563 ; */
 
    bool gmn = btp->get_attr_table().get_attr("grid_mapping_name") == "latitude_longitude";
    bool lpm = btp->get_attr_table().get_attr("longitude_of_prime_meridian") == "0.0";
    bool sma = btp->get_attr_table().get_attr("semi_major_axis") == "6378137.0";
    bool iflat = btp->get_attr_table().get_attr("inverse_flattening") == "298.257223563";
 
    return gmn && lpm && sma && iflat;
}

string FONgGrid::get_projection(DDS *dds)
{
    // Here's the information about the CF and projections
    // http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#grid-mappings-and-projections
    // How this code looks for mapping information: Look for an
    // attribute named 'grid_mapping' and get it's value. This attribute
    // with be the name of a variable in the dataset, so get that
    // variable. Now, look at the attributes of that variable.
    string mapping_info = d_grid->get_attr_table().get_attr("grid_mapping");
    if (mapping_info.empty()) mapping_info = d_grid->get_array()->get_attr_table().get_attr("grid_mapping");
 
    string WK_GCS = GeoTiffTransmitter::default_gcs;
 
    if (!mapping_info.empty()) {
        // "WGS84": same as "EPSG:4326" but has no dependence on EPSG data files.
        // "WGS72": same as "EPSG:4322" but has no dependence on EPSG data files.
        // "NAD27": same as "EPSG:4267" but has no dependence on EPSG data files.
        // "NAD83": same as "EPSG:4269" but has no dependence on EPSG data files.
        // "EPSG:n": same as doing an ImportFromEPSG(n).
 
        // The mapping info is actually stored as attributes of an Int32 variable.
        BaseType *btp = dds->var(mapping_info);
        if (btp && btp->name() == "crs") {
            if (is_wgs84(btp))
                WK_GCS = "WGS84";
            else if (is_spherical(btp)) WK_GCS = "EPSG:4047";
        }
    }
 
    OGRSpatialReference srs;
    srs.SetWellKnownGeogCS(WK_GCS.c_str());
    char *srs_wkt = NULL;
    srs.exportToWkt(&srs_wkt);
 
    string wkt = srs_wkt;   // save the result
 
    CPLFree(srs_wkt);       // free memory alloc'd by GDAL
 
    return wkt;
}
 
double *FONgGrid::get_data()
{
    if (!d_grid->get_array()->read_p()) d_grid->get_array()->read();
 
    // This code assumes read() has been called.
    return extract_double_array(d_grid->get_array());
}