ugr5.cc

Go to the documentation of this file.

// -*- mode: c++; c-basic-offset:4 -*-

// This file is part of libdap, A C++ implementation of the OPeNDAP Data
// Access Protocol.

// Copyright (c) 2002,2003,2011,2012 OPeNDAP, Inc.
// Authors: Nathan Potter <ndp@opendap.org>
//          James Gallagher <jgallagher@opendap.org>
//          Scott Moe <smeest1@gmail.com>
//          Bill Howe <billhowe@cs.washington.edu>
//
// 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 OPeNDAP, Inc. at PO Box 112, Saunderstown, RI. 02874-0112.

// NOTE: This file is built only when the gridfields library is linked with
// the netcdf_handler (i.e., the handler's build is configured using the
// --with-gridfields=... option to the 'configure' script).

#include "config.h"

#include <limits.h>

#include <cstdlib>      // used by strtod()
#include <cerrno>
#include <cmath>
#include <iostream>
#include <sstream>
//#include <cxxabi.h>

#include <curl/curl.h>

#define DODS_DEBUG

#include "BaseType.h"
#include "Int32.h"
#include "Str.h"
#include "Array.h"
#include "Structure.h"
#include "Error.h"

#include "BESDebug.h"
#include "util.h"

#include "ugrid_utils.h"
#include "MeshDataVariable.h"
#include "TwoDMeshTopology.h"
#include "NDimensionalArray.h"

#include "ugr5.h"


using namespace std;
using namespace libdap;

namespace ugrid {

static string ugrSyntax =
                "ugr5(dim:int32, rangeVariable:string, [rangeVariable:string, ... ] condition:string)";

struct UgridRestrictArgs {
        locationType dimension;
        vector<libdap::Array *> rangeVars;
        string filterExpression;
};


static void addRangeVar(DDS *dds, libdap::Array *rangeVar, map<string, vector<MeshDataVariable *> *> *rangeVariables) {

        MeshDataVariable *mdv = new MeshDataVariable();
        mdv->init(rangeVar);
        string meshVarName = mdv->getMeshName();

    BaseType *meshVar = dds->var(meshVarName);

    if(meshVar == 0){
        string msg = "The range variable '"+mdv->getName()+"' references the mesh variable '"+meshVarName+
                "' which cannot be located in this dataset.";
        BESDEBUG("ugrid", "addRangeVar() - " << msg  << endl);
        throw Error(no_such_variable,msg);
    }


        // Get the rangeVariable vector for this mesh name from the map.
        vector<MeshDataVariable *> *requestedRangeVarsForMesh;
        map<string, vector<MeshDataVariable *> *>::iterator mit = rangeVariables->find(meshVarName);
        if(mit == rangeVariables->end()){
                // Not there? Make a new one.
                BESDEBUG("ugrid", "addRangeVar() - MeshTopology object for '" << meshVarName <<"' does NOT exist. Getting a 'new' one... "  << endl);

                requestedRangeVarsForMesh =  new vector<MeshDataVariable *>();
                (*rangeVariables)[meshVarName] = requestedRangeVarsForMesh;
        }
        else {
                // Sweet! Found it....
                BESDEBUG("ugrid", "addRangeVar() - MeshTopology object for '" << meshVarName <<"' exists. Retrieving... "  << endl);
                requestedRangeVarsForMesh = mit->second;
        }

        requestedRangeVarsForMesh->push_back(mdv);

}

static UgridRestrictArgs processUgrArgs(int argc, BaseType *argv[]) {

        BESDEBUG( "ugrid", "processUgrArgs() - BEGIN" << endl);

        UgridRestrictArgs args;
        args.rangeVars = vector<libdap::Array *>();

        // Check number of arguments; BESDEBUG is a macro. Use #define
        // DODS_DEBUG to activate the debugging stuff.
        if (argc < 3)
                throw Error(malformed_expr,
                                "Wrong number of arguments to ugrid restrict function: "
                                                + ugrSyntax + " was passed "
                                                + long_to_string(argc) + " argument(s)");

        BaseType * bt;

        // ---------------------------------------------
        // Process the first arg, which is the rank of the Restriction Clause
        bt = argv[0];
        if (bt->type() != dods_int32_c)
                throw Error(malformed_expr,
                                "Wrong type for first argument, expected DAP Int32. "
                                                + ugrSyntax + "  was passed a/an "
                                                + bt->type_name());
        args.dimension = (locationType) dynamic_cast<Int32&>(*argv[0]).value();
    BESDEBUG( "ugrid", "args.dimension: " << libdap::long_to_string(args.dimension) << endl);


        // ---------------------------------------------
        // Process the last argument, the relational/filter expression used to restrict the ugrid content.
        bt = argv[argc - 1];
        if (bt->type() != dods_str_c)
                throw Error(malformed_expr,
                                "Wrong type for third argument, expected DAP String. "
                                                + ugrSyntax + "  was passed a/an "
                                                + bt->type_name());
        args.filterExpression = dynamic_cast<Str&>(*bt).value();
    BESDEBUG( "ugrid", "args.filterExpression: '" << args.filterExpression << "' (AS RECEIVED)" << endl);

        int decodedLength;
        CURL *curl = curl_easy_init( );
        char *decodedFilterExpression = curl_easy_unescape(curl , args.filterExpression.c_str() ,  args.filterExpression.size() , &decodedLength );
        curl_easy_cleanup(curl );
        args.filterExpression = string(decodedFilterExpression,decodedLength);
    BESDEBUG( "ugrid", "args.filterExpression: '" << args.filterExpression << "' (URL DECODED)" << endl);
        curl_free(decodedFilterExpression);

        // --------------------------------------------------
        // Process the range variables selected by the user.
        // We know that argc>=3, because we checked so the
        // following loop will try to find at least one rangeVar,
        // and it won't try to process the first or last members
        // of argv.
        for (int i = 1; i < (argc - 1); i++) {
                bt = argv[i];
                if (bt->type() != dods_array_c)
                        throw Error(malformed_expr,
                                        "Wrong type for second argument, expected DAP Array. "
                                                        + ugrSyntax + "  was passed a/an "
                                                        + bt->type_name());

                libdap::Array *newRangeVar = dynamic_cast<libdap::Array*>(bt);
                if (newRangeVar == 0) {
                        throw Error(malformed_expr,
                                        "Wrong type for second argument. " + ugrSyntax
                                                        + "  was passed a/an " + bt->type_name());
                }
                args.rangeVars.push_back(newRangeVar);
        }

        BESDEBUG("ugrid", "processUgrArgs() - END" << endl);

        return args;

}


static string arrayState(libdap::Array *dapArray, string indent){

    libdap::Array::Dim_iter thisDim;
    stringstream arrayState;
    arrayState << endl;
    arrayState << indent << "dapArray: '" << dapArray->name() << "'   " ;
    arrayState << "type: '" << dapArray->type_name( ) << "'  " ;

    arrayState << "read(): " << dapArray->read() << "  ";
    arrayState << "read_p(): " << dapArray->read_p() << "  ";
    arrayState << endl;

    for(thisDim= dapArray->dim_begin(); thisDim!=dapArray->dim_end() ;thisDim++){

    arrayState << indent << "  dim:  '" << dapArray->dimension_name(thisDim) << "' ";
    arrayState << indent << "    start:  " << dapArray->dimension_start(thisDim)  << " ";
    arrayState << indent << "    stride: " << dapArray->dimension_stride(thisDim)  << " ";
    arrayState << indent << "    stop:   " << dapArray->dimension_stop(thisDim)  << " ";
    arrayState << endl;
    }
    return arrayState.str();
}


static void copyUsingSubsetIndex(libdap::Array *sourceArray, vector<unsigned int> *subsetIndex, void *result){
    BESDEBUG("ugrid", "ugrid::copyUsingIndex() - BEGIN" << endl);

    switch (sourceArray->var()->type()) {
    case dods_byte_c:
        sourceArray->value(subsetIndex, (dods_byte *)result);
        break;
    case dods_uint16_c:
        sourceArray->value(subsetIndex, (dods_uint16 *)result);
        break;
    case dods_int16_c:
        sourceArray->value(subsetIndex, (dods_int16 *)result);
        break;
    case dods_uint32_c:
        sourceArray->value(subsetIndex, (dods_uint32 *)result);
        break;
    case dods_int32_c:
        sourceArray->value(subsetIndex, (dods_int32 *)result);
        break;

    case dods_float32_c:
        sourceArray->value(subsetIndex, (dods_float32 *)result);
        break;
    case dods_float64_c:
        sourceArray->value(subsetIndex, (dods_float64 *)result);
        break;

    default:
        throw InternalErr(__FILE__, __LINE__,
                "ugrid::hgr5::copyUsingSubsetIndex() - Unknown DAP type encountered.");
    }

    BESDEBUG("ugrid", "ugrid::copyUsingIndex() - END" << endl);
}

static string vectorToString(vector<unsigned int> *index){
    BESDEBUG("ugrid", "indexToString() - BEGIN"<< endl);
    BESDEBUG("ugrid", "indexToString() - index.size(): " << libdap::long_to_string(index->size()) << endl);
    stringstream s;
    s << "[";

    for(unsigned int i=0; i<index->size() ; ++i){
        s  << ((i>0)?", ":" ");
        s  << (*index)[i];
    }
    s << "]";
    BESDEBUG("ugrid", "indexToString() - END"<< endl);

    return s.str();
}




static void rDAWorker(MeshDataVariable *mdv, libdap::Array::Dim_iter thisDim, vector<unsigned int> *slab_subset_index,
                NDimensionalArray *results)
{
        libdap::Array *dapArray = mdv->getDapArray();

        BESDEBUG("ugrid", "rDAWorker() - slab_subset_index" << vectorToString(slab_subset_index) << " size: "
                << libdap::long_to_string(slab_subset_index->size()) << endl);

        // The locationCoordinateDimension is the dimension of the array that is associated with the ugrid "rank" - e.g. it is the
        // dimension that ties the variable to the 'nodes' (rank 0) or 'edges' (rank 1) or 'faces' (rank 2) of the ugrid.
        libdap::Array::Dim_iter locationCoordinateDim = mdv->getLocationCoordinateDimension();

        BESDEBUG("ugrid", "rdaWorker() - thisDim: '" << dapArray->dimension_name(thisDim) << "'" << endl);
        BESDEBUG("ugrid", "rdaWorker() - locationCoordinateDim: '" << dapArray->dimension_name(locationCoordinateDim) << "'" << endl);
        // locationCoordinateDim is, e.g., 'nodes'. jhrg 10/25/13
        if (thisDim != locationCoordinateDim) {
                unsigned int start = dapArray->dimension_start(thisDim, true);
                unsigned int stride = dapArray->dimension_stride(thisDim, true);
                unsigned int stop = dapArray->dimension_stop(thisDim, true);

                BESDEBUG("ugrid", "rdaWorker() - array state: " << arrayState(dapArray, "    "));

                for (unsigned int dimIndex = start; dimIndex <= stop; dimIndex += stride) {
                        dapArray->add_constraint(thisDim, dimIndex, 1, dimIndex);
                        rDAWorker(mdv, thisDim + 1, slab_subset_index, results);
                }

                // Reset the constraint for this dimension.
                dapArray->add_constraint(thisDim, start, stride, stop);
        }
        else {
                BESDEBUG("ugrid", "rdaWorker() - Found locationCoordinateDim" << endl);

                if ((thisDim + 1) != dapArray->dim_end()) {
                        string msg = "rDAWorker() - The location coordinate dimension is not the last dimension in the array. Hyperslab subsetting of this dimension is not supported.";
                        BESDEBUG("ugrid", msg << endl);
                        throw Error(malformed_expr, msg);
                }

                BESDEBUG("ugrid", "rdaWorker() - Arrived At Last Dimension" << endl);

                dapArray->set_read_p(false);

                BESDEBUG("ugrid", "rdaWorker() - dap array: " << arrayState(dapArray, "    "));

                vector<unsigned int> lastDimHyperSlabLocation;
                NDimensionalArray::retrieveLastDimHyperSlabLocationFromConstrainedArrray(dapArray, &lastDimHyperSlabLocation);

                BESDEBUG("ugrid", "rdaWorker() - lastDimHyperSlabLocation: "
                        << NDimensionalArray::vectorToIndices(&lastDimHyperSlabLocation) << endl);

                // unused. 4/7/14 jhrg. unsigned int elementCount;

                void *slab;
        //results->getLastDimensionHyperSlab(&lastDimHyperSlabLocation, &slab, &elementCount);
        results->getNextLastDimensionHyperSlab(&slab);

                dapArray->read();

                copyUsingSubsetIndex(dapArray, slab_subset_index, slab);
        }
}

static libdap::Array *restrictRangeVariableByOneDHyperSlab(
        MeshDataVariable *mdv,
        vector<unsigned int> *slab_subset_index
        ){

    long restrictedSlabSize = slab_subset_index->size();

    BESDEBUG("ugrid", "restrictRangeVariableByOneDHyperSlab() - slab_subset_index"<< vectorToString(slab_subset_index) <<
            " size: " <<  libdap::long_to_string(restrictedSlabSize) << endl);

    libdap::Array *sourceDapArray = mdv->getDapArray();

    BESDEBUG("ugrid", "restrictDapArrayByOneDHyperSlab() - locationCoordinateDim: '" <<
            sourceDapArray->dimension_name(mdv->getLocationCoordinateDimension()) << "'" << endl);

    // We want the manipulate the Array's Dimensions so that only a single dimensioned slab of the location coordinate dimension
    // is read at a time. We need to cache the original constrained dimensions so that we can build the correct collection of
    // location coordinate dimension slabs.

    // Now we need to compute the shape of the final subset result array from the source range variable array and the slab subset.
    // What's the shape of the source array with any constraint applied?
    vector<unsigned int> resultArrayShape(sourceDapArray->dimensions(true));
    NDimensionalArray::computeConstrainedShape(sourceDapArray, &resultArrayShape );

    stringstream msg;
    for(unsigned int i=0; i< resultArrayShape.size(); i++){
        msg << "[" << resultArrayShape[i] << "]";
    }
    BESDEBUG("ugrid", "restrictDapArrayByOneDHyperSlab() - Constrained source array shape" << msg.str() << endl);
    msg.str("");

    // Now, we know that the result array has a last dimension slab size determined by the slab_subset_index (whichwas made by the ugrid sub-setting),
    // so we make the result array shape reflect that

    resultArrayShape[sourceDapArray->dimensions(true) - 1] = restrictedSlabSize; // jhrg 10/25/13 resultArrayShape.last() = ...
    libdap::Type dapType = sourceDapArray->var()->type();


    BESDEBUG("ugrid", "restrictDapArrayByOneDHyperSlab() - UGrid restricted HyperSlab has  "<< restrictedSlabSize << " elements." << endl);
    BESDEBUG("ugrid", "restrictDapArrayByOneDHyperSlab() - Array is of type '"<< libdap::type_name(dapType) << "'" << endl);

    for(unsigned int i=0; i< resultArrayShape.size(); i++){
        msg << "[" << resultArrayShape[i] << "]";
    }
    BESDEBUG("ugrid", "restrictDapArrayByOneDHyperSlab() - resultArrayShape" << msg.str() << endl);
    msg.str(std::string());

    // Now we make a new NDimensionalArray instance that we will use to hold the results.
    NDimensionalArray *result = new NDimensionalArray(&resultArrayShape, dapType);

    // And we pass that along with other stuff into the recursive rDAWorker that's going to go get all the stuff
    rDAWorker(mdv, sourceDapArray->dim_begin(), slab_subset_index, result);

    // And now that the recursion we grab have the NDimensionalArray cough up the rteuslt as a libdap::Array
    libdap::Array  *resultDapArray = result->getArray(sourceDapArray);

    // Delete the NDimensionalArray
    delete result;

    // Return the result as libdap:Array
    return resultDapArray;
}

void ugr5(int argc, BaseType *argv[], DDS &dds, BaseType **btpp)
{
        BESDEBUG("ugrid", "ugr5() - BEGIN" << endl);

        string info = string("<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n")
                        + "<function name=\"ugr5\" version=\"0.1\">\n"
                        + "Server function for Unstructured grid operations.\n" + "usage: "
                        + ugrSyntax + "\n"+ "</function>";

        if (argc == 0) {
                Str *response = new Str("info");
                response->set_value(info);
                *btpp = response;
                return;
        }

        // Process and QC the arguments
        UgridRestrictArgs args = processUgrArgs(argc, argv);

        // Each range variable is associated with a "mesh" i.e. a mesh topology variable. Since there may be more than one mesh in a
        // dataset, and the user may request more than one range variable for each mesh we need to sift through the list of requested
        // range variables and organize them by mesh topology variable name.
        map<string, vector<MeshDataVariable *> *> *meshToRangeVarsMap = new map<string, vector<MeshDataVariable *> *>();

        // For every Range variable in the arguments list, locate it and ingest it.
        vector<libdap::Array *>::iterator it;
        for (it = args.rangeVars.begin(); it != args.rangeVars.end(); ++it) {
            addRangeVar(&dds, *it, meshToRangeVarsMap);
        }
        BESDEBUG("ugrid", "ugr5() - The user requested "<< args.rangeVars.size() << " range data variables." << endl);
        BESDEBUG("ugrid", "ugr5() - The user's request referenced "<< meshToRangeVarsMap->size() << " mesh topology variables." << endl);

        // ----------------------------------
        // OK, so up to this point we have not read any data from the data set, but we have QC'd the inputs and verified that
        // it looks like the request is consistent with the semantics of the dataset.
        // Now it's time to read some data and pack it into the GridFields library...


    // TODO This returns a single structure but it would make better sense to the
    // world if it could return a vector of objects and have them appear at the
    // top level of the DDS.
    // FIXME fix the names of the variables in the mesh_topology attributes
    // If the server side function can be made to return a DDS or a collection of BaseType's then the
    // names won't change and the original mesh_topology variable and it's metadata will be valid
        Structure *dapResult = new Structure("ugr_result");


        // Since we only want each ugrid structure to appear in the results one time  (cause otherwise we might be trying to add
        // the same variables with the same names to the result multiple times.) we grind on this by iterating over the
        // names of the mesh topology names.
    vector<MeshDataVariable *>::iterator rvit;
    map<string, vector<MeshDataVariable *> *>::iterator mit;
        for (mit = meshToRangeVarsMap->begin(); mit != meshToRangeVarsMap->end(); ++mit) {

                string meshVariableName = mit->first;
                vector<MeshDataVariable *> *requestedRangeVarsForMesh = mit->second;

        vector<BaseType *> dapResults;

        // Building the restricted TwoDMeshTopology without adding any range variables and then converting the result
        // Grid field to Dap Objects should return all of the Ugrid structural stuff - mesh variable, node coordinate variables,
        // face and edge coordinate variables if present.
        BESDEBUG("ugrid", "ugr5() - Adding restricted mesh_topology structure for mesh '" << meshVariableName << "' to DAP response." << endl);

        TwoDMeshTopology *tdmt = new TwoDMeshTopology();
        tdmt->init(meshVariableName, &dds);

        tdmt->buildBasicGfTopology();
        tdmt->addIndexVariable(node);
        tdmt->addIndexVariable(face);
        tdmt->applyRestrictOperator(args.dimension, args.filterExpression);

        // 3: because there are nodes (rank = 0), edges (rank = 1), and faces (rank = 2). jhrg 10/25/13
        vector<vector<unsigned int> *> location_subset_indices(3);

        long nodeResultSize = tdmt->getResultGridSize(node);
        BESDEBUG("ugrid", "ugr5() - there are "<< nodeResultSize << " nodes in the subset." << endl);
        vector<unsigned int> node_subset_index(nodeResultSize);
        if(nodeResultSize>0){
                tdmt->getResultIndex(node, &node_subset_index[0]);
        }
        location_subset_indices[node] = &node_subset_index;

        BESDEBUG("ugrid", "ugr5() - node_subset_index"<< vectorToString(&node_subset_index) << endl);


        long faceResultSize = tdmt->getResultGridSize(face);
        BESDEBUG("ugrid", "ugr5() - there are "<< faceResultSize << " faces in the subset." << endl);
        vector<unsigned int> face_subset_index(faceResultSize);
        if(faceResultSize > 0){
                tdmt->getResultIndex(face, &face_subset_index[0]);
        }
        location_subset_indices[face] = &face_subset_index;
        BESDEBUG("ugrid", "ugr5() - face_subset_index: "<< vectorToString(&face_subset_index) << endl);



        // This gets all the stuff that's attached to the grid - which at this point does not include the range variables but does include the
        // index variable. good enough for now but need to drop the index....
        tdmt->convertResultGridFieldStructureToDapObjects(&dapResults);


        BESDEBUG("ugrid", "ugr5() - Restriction of mesh_topology '"<< tdmt->getMeshVariable()->name() << "' structure completed." << endl);

        // now that we have the mesh topology variable we are going to look at each of the requested
        // range variables (aka MeshDataVariable instances) and we're going to subset that using the
        // gridfields library and add its subset version to the results.
            for(rvit=requestedRangeVarsForMesh->begin(); rvit!=requestedRangeVarsForMesh->end(); rvit++){
                MeshDataVariable *mdv = *rvit;

                BESDEBUG("ugrid", "ugr5() - Processing MeshDataVariable  '"<< mdv->getName() << "' associated with rank/location: "<< mdv->getGridLocation() << endl);


                tdmt->setLocationCoordinateDimension(mdv);



            libdap::Array *restrictedRangeVarArray =
                    restrictRangeVariableByOneDHyperSlab(mdv, location_subset_indices[mdv->getGridLocation()]);

            BESDEBUG("ugrid", "ugr5() - Adding resulting dapArray  '"<< restrictedRangeVarArray->name() << "' to dapResults." << endl);

            dapResults.push_back(restrictedRangeVarArray);


            }

        delete tdmt;


                BESDEBUG("ugrid", "ugr5() - Adding GF::GridField results to DAP structure " << dapResult->name() << endl);
                for (vector<BaseType *>::iterator btIt=dapResults.begin(); btIt != dapResults.end(); ++btIt) {
                    BaseType *bt = *btIt;
                BESDEBUG("ugrid", "ugr5() - Adding variable "<< bt->name() << " to DAP structure " << dapResult->name() << endl);
                        dapResult->add_var_nocopy(bt);
                }
        }

        *btpp = dapResult;

        BESDEBUG("ugrid", "ugr5() - Releasing maps and vectors..." << endl);
        for (mit = meshToRangeVarsMap->begin(); mit != meshToRangeVarsMap->end(); ++mit) {
            vector<MeshDataVariable *> *requestedRangeVarsForMesh = mit->second;
        for(rvit=requestedRangeVarsForMesh->begin(); rvit!=requestedRangeVarsForMesh->end(); rvit++){
            MeshDataVariable *mdv = *rvit;
            delete mdv;
        }
                delete requestedRangeVarsForMesh;
        }
    delete meshToRangeVarsMap;






        BESDEBUG("ugrid", "ugr5() - END" << endl);

        return;
}


} // namespace ugrid_restrict