Logo Search packages:      
Sourcecode: magics++ version File versions  Download package

GribRegularInterpretor.cc

Go to the documentation of this file.
/******************************** LICENSE ********************************

 Copyright 2007 European Centre for Medium-Range Weather Forecasts (ECMWF)

 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at 

    http://www.apache.org/licenses/LICENSE-2.0

 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.

 ******************************** LICENSE ********************************/

/*! \file GribRegularInterpretor.cc
    \brief Implementation of the Template class GribRegularInterpretor.
    \author Graphics Section, ECMWF

    Started: Mon 18-Apr-2005

    Changes:
*/

#include <limits>

#include "TeProjection.h"
#include "GribRegularInterpretor.h"
#include "GribDecoder.h"
#include "Matrix.h"
#include "LocalTable.h"
#include "RasterData.h"


using namespace magics;

GribRegularInterpretor::GribRegularInterpretor() 
{
}
GribRegularInterpretor::~GribRegularInterpretor() 
{
}


void GribInterpretor::scaling(const GribDecoder& grib, Matrix** matrix) const
{

// Fisrt check that they are not revide fields! 

      long derived = grib.getLong("generatingProcessIdentifier");
      
      
      
      
      
      // later... 
      long param = grib.getLong("indicatorOfParameter");
      long centre = grib.getLong("identificationOfOriginatingGeneratingCentre");
      long table = grib.getLong("gribTablesVersionNo");
      
      if (grib.getScaling())
      {
            if (derived != 254)  {
            try {
                  const ParamDef& paramdef = LocalTable::localInfo(param, table, centre);
                  (*matrix)->multiply(paramdef.scaling());
                  (*matrix)->plus(paramdef.offset());
            }
            catch (...)
            {
                  Log::warning() << " Can not find information for the parameter [" << param << "." << table << "]\n";
            }
            }
      }
      else {
            (*matrix)->multiply(grib.getScaling_factor());
            (*matrix)->plus(grib.getScaling_offset());
      }

}

/*!
 Class information are given to the output-stream.
*/          
void GribRegularInterpretor::print(ostream& out)  const
{
      out << "GribRegularInterpretor[";
      out << "]";
}


void GribRegularInterpretor::interpretAsMatrix(const GribDecoder& grib, Matrix** matrix) const 
{


  Log::dev() << "GribRegularInterpretor::interpretAsMatrix" << "\n";
  long nblon = grib.getLong("numberOfPointsAlongAParallel");
  long nblat = grib.getLong("numberOfPointsAlongAMeridian");
  
  if ( *matrix == 0 ) *matrix = new Matrix(nblat, nblon);
  
  
  size_t nb;
  grib_get_size(grib.id(), "values", &nb);

  

  Log::dev() << "numberOfFieldValues[" << nb << "]" << "\n";
  double missing = -std::numeric_limits<double>::max();
  grib.setDouble("missingValue", missing);    
  (*matrix)->missing(missing);
  (*matrix)->akimaEnabled();

  double north = grib.getDouble("latitudeOfFirstGridPointInDegrees");
  double west = grib.getDouble("longitudeOfFirstGridPointInDegrees");
  double south = grib.getDouble("latitudeOfLastGridPointInDegrees");; 
  double east = grib.getDouble("longitudeOfLastGridPointInDegrees");;
  

  Log::dev() << "NewAPI---> area[" << west << ", " << north << ", " << east << ", " << south << "]" << "\n";
  double loni = longitudeIncrement(grib);
  
  double lon = (east-west)/(nblon-1);
  
  Log::dev() << "increment -->" << loni << " (from->" << west << " to-->" << west + (nblon-1) *loni << ")" <<  endl;
  Log::dev() << "calcul -->" << lon << " (from->" << west << " to-->" << west + (nblon-1) *lon << ")" <<  endl;

  latitudes(grib, (*matrix)->rowsAxis());
  
  
 
 

 
   double x = west;
      for (int i = 0; i < nblon; i++)
      {
          
            (*matrix)->columnsAxis().push_back(x);
            x  = west + (i+1)*lon;
            
      }
      
 
  
  (*matrix)->setMapsAxis();

  
 try { 
       (*matrix)->resize(nb);
      size_t aux = size_t(nb);
      grib_get_double_array(grib.id(),"values", &(*matrix)->front(),&aux); 
      (*matrix)->missing(missing);
 } 
 catch (...) 
 {
       throw MagicsException("Not enough memory");
 }
  
}

void  GribRegularInterpretor::latitudes(const GribDecoder& grib, vector<double>& latitudes) const
{
      
            double north = grib.getDouble("latitudeOfFirstGridPointInDegrees");     

            long nblat = grib.getLong("numberOfPointsAlongAMeridian");
            int scanning = grib.getLong("jScansPositively") ? 1 : -1;
            double lat =  scanning * grib.getDouble("jDirectionIncrementInDegrees");
      
            double y = north;
            for (int i = 0; i < nblat; i++)
            {
                     
                        latitudes.push_back(y);
                         y  += lat;
            }

}

void GribRegularGaussianInterpretor::latitudes(const GribDecoder& grib, vector<double>& latitudes) const 
{
      long res     = grib.getLong("numberOfParallelsBetweenAPoleAndTheEquator");
      double array[2 *res];
      grib_get_gaussian_latitudes(res, array);

      for ( int i = 0; i < 2*res; i++ )
      {
                  latitudes.push_back(array[i]);
      }
      
}

double GribRegularInterpretor::longitudeIncrement(const GribDecoder& grib) const 
{
      int scanning = grib.getLong("iScanNegatively") ? -1 : 1;
      return scanning * grib.getDouble("iDirectionIncrementInDegrees"); 
      /*
      double north = grib.getDouble("latitudeOfFirstGridPointInDegrees");     
      double south = grib.getDouble("latitudeOfLastGridPointInDegrees");; 

      long nblat = grib.getLong("numberOfPointsAlongAMeridian");

      // here we have to calculate the latitude increment 
      return (north - south)/(nblat-1);
      */
}




void GribRegularInterpretor::interpretAsRaster(const GribDecoder& grib, RasterData<GeoPoint>& raster) const
{
      Log::dev() << "GribRegularInterpretor::interpretAsRaster" << "\n";
      
      
      
       double north = grib.getDouble("latitudeOfFirstGridPointInDegrees");
        double west = grib.getDouble("longitudeOfFirstGridPointInDegrees");
        double south = grib.getDouble("latitudeOfLastGridPointInDegrees");; 
        double east = grib.getDouble("longitudeOfLastGridPointInDegrees");;
      Log::dev() << "area[" << north << ", " << west << ", " << south << ", " << east << "]" << "\n";
      
      raster.setUpperRightCorner(north, east);
      raster.setLowerLeftCorner(south, west);

      
      double lon = longitudeIncrement(grib);
      int scanning = grib.getLong("jScansPositively") ? 1 : -1;
      double lat =  scanning * grib.getDouble("jDirectionIncrementInDegrees");
      
      
      
      raster.setXResolution(lon);
      raster.setYResolution(lat);
      
       long nblon = grib.getLong("numberOfPointsAlongAParallel");
        long nblat = grib.getLong("numberOfPointsAlongAMeridian");
      
      raster.setColumns(nblon);
      raster.setRows(nblat);
//    
      raster.setProjection(new TeLatLong(TeDatum()));
      double missing = std::numeric_limits<double>::max();
      size_t nb;
      grib_get_size(grib.id(), "values", &nb);
 
      double *data = new double[nb];
        
      size_t aux = size_t(nb);
      
      grib_get_double_array(grib.id(),"values",data,&aux); 
      vector<double> tmp;
       for (int i = 0; i < nb; i++)
              tmp.push_back(data[i]);
       
       double min = *(std::min_element(tmp.begin(), tmp.end()));
       double max = *(std::max_element(tmp.begin(), tmp.end()));
  
       double step = (max-min)/126;
       
    raster.reserve(nb);
      
    for (int i = 0; i < nb; i++) {
      
      //cout << data[i] << "---> " << int ((data[i]/step) ) << endl;
        raster.push_back(data[i]);
    } 
}


void GribReducedGaussianInterpretor::print(ostream& out)  const
{
      out << "GribRegularInterpretor[";
      out << "]";
}

void GribReducedGaussianInterpretor::interpretAsMatrix(const GribDecoder& grib, Matrix** matrix) const
{
      Log::dev() << "GribRegularInterpretor::interpretAsMatrix" << "\n";
    Log::dev() << "GribRegularInterpretor::interpretAsMatrix" << "\n";
 
  
    *matrix = new Matrix();
      size_t nb;
      grib_get_size(grib.id(), "values", &nb);

      Log::dev() << "numberOfFieldValues[" << nb << "]" << "\n";
      double missing = -std::numeric_limits<double>::max();
      grib.setDouble("missingValue", missing);
      (*matrix)->missing(missing);
      (*matrix)->akimaEnabled();

      double north = grib.getDouble("latitudeOfFirstGridPointInDegrees");
      double west  = grib.getDouble("longitudeOfFirstGridPointInDegrees");
      double south = grib.getDouble("latitudeOfLastGridPointInDegrees");; 
      double east  = grib.getDouble("longitudeOfLastGridPointInDegrees");;
      double plp   = grib.getDouble("PLPresent");
      long res     = grib.getLong("numberOfParallelsBetweenAPoleAndTheEquator");

      Log::dev() << "NewAPI---> area[" << west << ", " << north << ", " << east << ", " << south << "]" << "\n";
      Log::dev() << "PLPresent---> " << plp << "\n";
      Log::dev() << "Res---> " << res << "\n";

      double pl[2*res];

      size_t aux = 2*res;
      grib_get_double_array(grib.id(),"pl",pl,&aux);
      // We have to determine if the field is global! 
      bool global = same(north-south, 180.);
      if (global) Log::dev() << "YES THE FIELD IS GLOBAL" << endl;

      double *data = new double[nb];

      size_t aux2 = size_t(nb);
      int nblon = 4*res;
      double width = east-west;
      double step = (width)/(nblon-1); 

      grib_get_double_array(grib.id(),"values",data,&aux2);
      int d = 0;
      for ( size_t i = 0; i < aux; i++)
      {
            vector<double> p;
            for ( int ii = 0; ii < pl[i]; ii++)
            {
                  p.push_back(data[d]);
                  d++;
            }
            
            if (global) 
                  p.push_back(p.front());
            
            double lon = west;
            int p1 = 0;
            int p2 = 1;
            double lon1 = west;
            double lon2 = lon1 + (width/(p.size()-1));

            for ( int x = 0; x < 4*res; x++ ) {
                  double d1 = (lon2 - lon)/(lon2-lon1);
                  double d2 = 1-d1;
                  double val;
                  if (p[p1] == missing ||  p[p2] == missing)
                                    val = missing;
                              else
                                    val = (p[p1] * d1) + (p[p2] * d2);
                  (*matrix)->push_back(val);
                  lon += step;
                  if ( lon >= lon2 ) {
                        p1++;
                        p2++;
                        lon1 = lon2;
                        lon2 += (width)/(p.size()-1);
                  }
            }
      }
    
      delete [] data;


      
      double lon = (width) / (nblon-1);
      
            for (int x = 0; x < nblon; x++)
            {
                  (*matrix)->columnsAxis().push_back(west+(x*lon));
            }
            
      double array[2 *res];
      long par = grib.getLong("numberOfParallelsBetweenAPoleAndTheEquator");
      grib_get_gaussian_latitudes(par, array);

      for ( int i = 0; i < 2*res; i++ )
      {
            (*matrix)->rowsAxis().push_back(array[i]);
      }
      (*matrix)->setMapsAxis();
}

void GribReducedLatLonInterpretor::print(ostream& out)  const
{
      out << "GribReducedLatLonInterpretor[";
      out << "]";
}

void GribReducedLatLonInterpretor::interpretAsMatrix(const GribDecoder& grib, Matrix** matrix) const
{
      
 
  
    *matrix = new Matrix();
      size_t nb;
      grib_get_size(grib.id(), "values", &nb);

      Log::dev() << "numberOfFieldValues[" << nb << "]" << "\n";
      double missing = std::numeric_limits<double>::max();
      grib.setDouble("missingValue", missing);
      (*matrix)->missing(missing);
      (*matrix)->akimaEnabled();

      double north = grib.getDouble("latitudeOfFirstGridPointInDegrees");
      double west  = grib.getDouble("longitudeOfFirstGridPointInDegrees");
      double south = grib.getDouble("latitudeOfLastGridPointInDegrees");; 
      double east  = grib.getDouble("longitudeOfLastGridPointInDegrees");;
      
      size_t res     = grib.getLong("Nj");

      Log::dev() << "NewAPI---> area[" << west << ", " << north << ", " << east << ", " << south << "]" << "\n";
      Log::dev() << "Res---> " << res << "\n";

      double pl[res];
      
      long nblat = grib.getLong("numberOfPointsAlongAMeridian");
      int scanning = grib.getLong("jScansPositively") ? 1 : -1;
      double lat =  scanning * grib.getDouble("jDirectionIncrementInDegrees");

      grib_get_double_array(grib.id(),"pl",pl,&res);

      double *data = new double[nb];

      size_t aux2 = size_t(nb);
      int nblon = 2*res;
      float width = east-west;
      float step = (width)/(nblon-1); 
      // We have to determine if the field is global! 
      bool global = (same(east + abs(lat) - width, 0) ||  width + abs(lat) >= 360);
      if (global) Log::dev() << "YES THE FIELD IS GLOBAL" << endl;
      
      grib_get_double_array(grib.id(),"values",data,&aux2);
      int d = 0;
      for ( size_t i = 0; i < res; i++)
      {
            
            
            float lon = west;
      
            float lon1 = west;
            if ( pl[i] == 0 ) {
                  // add missing data 
                  for ( int x = 0; x < nblon; x++ ) 
                        (*matrix)->push_back(missing);
            }
            else {
                  int p1 = 0;
                  int p2 = 1;
                  vector<double> p;
                  for ( int ii = 0; ii < pl[i]; ii++)
                  {
                                    p.push_back(data[d]);
                                    d++;
                  }
            if (global) p.push_back(p.front());
            
            float lon2 = lon1 + (width/(p.size()-1));

            for ( int x = 0; x < nblon; x++ ) {
                  float d1 = (lon2 - lon)/(lon2-lon1);
                  float d2 = 1-d1;
                  double val; 
                  if (p[p1] == missing ||  p[p2] == missing)
                        val = missing;
                  else
                        val = (p[p1] * d1) + (p[p2] * d2);
                  (*matrix)->push_back(val);
                  lon += step;
                  if ( lon >= lon2 ) {
                        p1++;
                        p2++;
                        lon1 = lon2;
                        lon2 += (width)/(p.size()-1);
                  }
            }
            }
      }
    
      delete [] data;


      float lon = (width) / (nblon-1);
//    float lat = (north - south) / (2 *res);
      for (int x = 0; x < nblon; x++)
      {
            (*matrix)->columnsAxis().push_back(west+(x*lon));
      }
      
      

                  
            
                  double y = north;
                  for (long i = 0; i < nblat; i++)
                  {
                        (*matrix)->rowsAxis().push_back(y);
                              
                               y  += lat;
                  }

      (*matrix)->setMapsAxis();
}


GeoPoint GribRotatedInterpretor::unrotate(double lat, double lon) const
{
        const double cToRadian         = 3.14/180.0;
        double zradi  = 1./cToRadian;
        double zsycen = sin(cToRadian*(southPoleLat_+90.));
        double zcycen = cos(cToRadian*(southPoleLat_+90.));

        double zsxrot = sin(cToRadian*lon);
        double zcxrot = cos(cToRadian*lon);
        double zsyrot = sin(cToRadian*lat);
        double zcyrot = cos(cToRadian*lat);
        double zsyreg = zcycen*zsyrot + zsycen*zcyrot*zcxrot;
        zsyreg = max( min(zsyreg, +1.0), -1.0 );

        double pyreg = asin(zsyreg)*zradi;

        double zcyreg = cos(pyreg*cToRadian);
        double zcxmxc = (zcycen*zcyrot*zcxrot - zsycen*zsyrot)/zcyreg;
        zcxmxc = max( min(zcxmxc, +1.0), -1.0 );
        double zsxmxc = zcyrot*zsxrot/zcyreg;
        double zxmxc  = acos(zcxmxc)*zradi;
        if( zsxmxc < 0.0)
          zxmxc = -zxmxc;

        double pxreg = zxmxc + southPoleLon_;
       
        return GeoPoint( pxreg, pyreg );
}






void GribRotatedInterpretor::print(ostream& out) const
{
      out << "GribRotatedInterpretor[]";
}


GeoPoint GribLambertAzimutalInterpretor::unrotate(double lat, double lon) const
{
        
        return GeoPoint(lon, lat);
}




void GribLambertAzimutalInterpretor::interpretAsMatrix(const GribDecoder& grib, Matrix** matrix) const
{
  
        long nblon = grib.getLong("numberOfPointsAlongXAxis");
        long nblat = grib.getLong("numberOfPointsAlongYAxis");
        
        ProjectedMatrix* lambert = new ProjectedMatrix(nblat, nblon);
        
        if ( *matrix == 0 ) *matrix = lambert;
        
        size_t nb;
        grib_get_size(grib.id(), "values", &nb);

        

        Log::dev() << "numberOfFieldValues[" << nb << "]" << "\n";
        double missing = -std::numeric_limits<double>::max();
        missing = grib.getDouble("missingValue");    
        lambert->missing(missing);

        double north = grib.getDouble("latitudeOfFirstGridPointInDegrees");
        double west = grib.getDouble("longitudeOfFirstGridPointInDegrees");
        
        

        Log::dev() << "NewAPI---> area[" << west << ", " << north << "]" << "\n";
        
        //lambert->standardParallel(grib.getDouble("standardParallelInMicrodegrees")/1000000);
       //lambert->centralLongitude(grib.getDouble("centralLongitudeInMicrodegrees")/1000000);
        
        
        double lon = grib.getDouble("xDirectionGridLengthInMillimetres")/1000;
        double lat =   grib.getDouble("xDirectionGridLengthInMillimetres")/1000;
       // double x = lambert->x(west, north);
        //double y = lambert->y(west, north);
        //double newlat = lambert->real_row(y,x);
        //double newlon = lambert->real_column(y,x);
        Log::dev() << "New API--->GridIncrements[" << lon << ", " << lat << "]" << "\n"; 
        Log::dev() << "New API--->GridDimensions[" << nblon << ", " << nblat << "]" << "\n";
        //Log::dev() << "New API--->Grid[<" << west << ", " << north << "], [" << x << ", " << y << "], [" << newlon << ", " << newlat << "\n";
         

        
      

        
        double *data = new double[nb]; 
        double *latitudes = new double[nb]; 
        double *longitudes = new double[nb];
        size_t aux = size_t(nb);
        grib_get_double_array(grib.id(),"values",data,&aux); 
        grib_get_double_array(grib.id(),"latitudes",latitudes,&aux); 
        grib_get_double_array(grib.id(),"longitudes", longitudes,&aux); 
        lambert->missing(missing);
        
        for (unsigned int i = 0; i < nb; i++)
        {
              lambert->push_back(data[i]);
              if (longitudes[i] > 180) {
                  longitudes[i] -= 360;
              }
              lambert->rowsAxis().push_back(latitudes[i]);   
              lambert->columnsAxis().push_back(longitudes[i]);
        }
        
        Log::dev() << "top-left[" << latitudes[0] << ", " << longitudes[0] << "]" << endl;
        Log::dev() << "top-right[" << latitudes[nblon-1] << ", " <<  longitudes[nblon-1] << "]" << endl; 
        Log::dev() << "top-right[" << latitudes[nb-1] << ", " <<  longitudes[nb-1] << "]" << endl;


        delete [] data;  
        delete [] latitudes;
        delete [] longitudes;  
        
        lambert->setMapsAxis();
        Log::dev() << "top-left[" << lambert->real_row(0,0) << ", " << lambert->real_column(0,0) << "]" << endl;
        Log::dev() << "top-right[" << lambert->real_row(0,nblon-1) << ", " << lambert->real_column(0,nblon-1) << "]" << endl;
        Log::dev() << "bottom-right[" << lambert->real_row(nblat-1, nblon-1) << ", " << lambert->real_column(nblat-1,nblon-1) << "]" << endl;
        Log::dev() << "bottom-left[" << lambert->real_row(nblat-1, 0) << ", " << lambert->real_column(nblat-1, 0) << "]" << endl;

        Log::dev() << **matrix<< endl;
    
    
}



void GribLambertAzimutalInterpretor::print(ostream& out) const
{
      out << "GribLambertAzimutalInterpretor[]";
}

void GribRotatedInterpretor::interpretAsMatrix(const GribDecoder& grib, Matrix** matrix) const
{
    long nblon = grib.getLong("numberOfPointsAlongAParallel");
    long nblat = grib.getLong("numberOfPointsAlongAMeridian");
  
    *matrix = new ProjectedMatrix(nblat, nblon);
    
    
    ProjectedMatrix* helper = static_cast<ProjectedMatrix*>(*matrix);
   
      
      
      southPoleLat_  = grib.getDouble("latitudeOfSouthernPoleInDegrees"); 
      southPoleLon_  = grib.getDouble("longitudeOfSouthernPoleInDegrees");
        double missing = -std::numeric_limits<double>::max();
              missing = grib.getDouble("missingValue");    
    
       helper->missing(missing);
        if ( *matrix == 0 ) *matrix = helper;
    

      int error;
      grib_iterator*  iterator = grib_iterator_new(grib.handle(), 0, &error);
      if ( error) {
            Log::warning() << "Pb readind data" << endl;
            return;
      }
      
      double lat, lon, value;
      int i = 0; 
      while ( grib_iterator_next(iterator, &lat, &lon, &value) ) { 
            GeoPoint point = unrotate(lat,lon);
            
            helper->rowsArray().push_back(point.y());   
            helper->columnsArray().push_back(point.x());
            (*matrix)->push_back( value);
            i++;
      }
      
      helper->setMapsAxis();


}

Generated by  Doxygen 1.6.0   Back to index