Changeset 782

Timestamp:

11/12/15 16:33:03 (7 years ago)

Author:

mhnguyen

Message:

Reading attributes of curvilinear grid from file

+) Correct some minor bugs detecting type of grid
+) Use constant string for attributes conforming to CF convention
+) Add part of code to read attributes of curvilinear grid

Test
+) On Curie
+) test_remap passes

Location:

XIOS/trunk

Files:

• inputs/REMAP/iodef.xml (modified) (5 diffs)
• src/io/inetcdf4.cpp (modified) (16 diffs)
• src/io/inetcdf4.hpp (modified) (1 diff)
• src/io/inetcdf4_decl.cpp (added)
• src/io/inetcdf4_impl.hpp (added)
• src/io/nc4_data_input.cpp (modified) (7 diffs)
• src/io/nc4_data_input.hpp (modified) (2 diffs)
• src/io/netCdfInterface.cpp (modified) (2 diffs)
• src/io/netCdf_cf_constant.hpp (added)
• src/node/domain.cpp (modified) (2 diffs)
• src/test/test_new_features.f90 (modified) (6 diffs)

XIOS/trunk/inputs/REMAP/iodef.xml

@@ -18 +18 @@
    <file_definition type="one_file" par_access="collective" output_freq="1ts" output_level="10" enabled=".TRUE.">
      <file id="output" name="output">
-        <field field_ref="src_field" name="field" />
+<!--        <field field_ref="src_field" name="field" />-->
      </file>
      <file id="output_dst" name="output_dst" type="one_file">
@@ -24 +24 @@
      </file>
      <file id="out_dst_regular_pole" name="out_dst_regular_pole" type="one_file">
-        <field field_ref="dst_field_regular_pole" name="field" />
+<!--        <field field_ref="dst_field_regular_pole" name="field" />-->
      </file>
      <file id="output_dst_regular" name="output_dst_regular" type="one_file">
-        <field field_ref="dst_field_regular" name="field" />
+<!--        <field field_ref="dst_field_regular" name="field" />-->
      </file>
      <file id="output_dst_curvilinear" name="output_dst_curvilinear" type="one_file">
-        <field field_ref="tmp_field_1" operation="instant"/>
+<!--        <field field_ref="tmp_field_1" operation="instant"/>-->
      </file>
 
@@ -44 +44 @@
         <field id="src_field_curvilinear" name="field_A" grid_ref="src_grid_curvilinear" operation="instant"/>
      </file>
+     <file id="output_surface_curvi" name="output_surface_curvi" mode="read" type="one_file">
+        <field id="field_surface_curvi" name="field_A" domain_ref="domain_surface_curvi" axis_ref="axis_surface_curvi" operation="instant"/>
+     </file>
    </file_definition>
 
@@ -49 +52 @@
    <axis_definition>
      <axis id="src_axis_curvilinear" />
+   </axis_definition>
+   <axis_definition>
+     <axis id="axis_surface_curvi" />
    </axis_definition>
 
@@ -72 +78 @@
      </domain>
      <domain id="src_domain_curvilinear" type="rectilinear">
+       <generate_rectilinear_domain />
+     </domain>
+     <domain id="domain_surface_curvi" type="curvilinear">
        <generate_rectilinear_domain />
      </domain>

XIOS/trunk/src/io/inetcdf4.cpp

@@ -1 +1 @@
 #include "inetcdf4.hpp"
 #include "netCdfInterface.hpp"
+#include "netCdf_cf_constant.hpp"
 
 #include <boost/algorithm/string.hpp>
@@ -302 +303 @@
     {
       const StdString& attname = *it;
-      if (attname.compare(name) == 0)
+      if (attname.compare(0, name.size(), name) == 0)
         return retvalue;
       retvalue++;
@@ -326 +327 @@
     {
       const StdString& attname = *it;
-      if (attname.compare(name) == 0) return true;
+      if (attname.compare(0, name.size(), name) == 0) return true;
     }
     return false;
@@ -339 +340 @@
     {
       const StdString& varname = *it;
-      if (varname.compare(name) == 0) return true;
+      if (varname.compare(0, name.size(), name) == 0) return true;
     }
     return false;
@@ -347 +348 @@
                                  const CVarPath* const path)
   {
-    return this->hasAttribute("coordinates", &name, path);
+    return this->hasAttribute(CCFKeywords::XIOS_CF_coordinates, &name, path);
   }
 
@@ -353 +354 @@
                             const CVarPath* const path)
   {
-    return this->hasAttribute("bounds", &name, path);
+    return this->hasAttribute(CCFKeywords::XIOS_CF_bounds, &name, path);
   }
 
@@ -441 +442 @@
   {
     StdString retvalue;
-    if (this->hasAttribute("coordinates", &name, path))
-    {
-      return this->getAttributeValue("coordinates", &name, path);
+    if (this->hasAttribute(CCFKeywords::XIOS_CF_coordinates, &name, path))
+    {
+      return this->getAttributeValue(CCFKeywords::XIOS_CF_coordinates, &name, path);
     }
     else
@@ -464 +465 @@
   {
     StdString retvalue;
-    if (this->hasAttribute("bounds", &name, path))
-      retvalue = this->getAttributeValue("bounds", &name, path);
+    if (this->hasAttribute(CCFKeywords::XIOS_CF_bounds, &name, path))
+      retvalue = this->getAttributeValue(CCFKeywords::XIOS_CF_bounds, &name, path);
     return retvalue;
   }
@@ -509 +510 @@
       return false;
 
+    bool isCurVi = true;
+    unsigned int nbLonLat = 0;
     std::list<StdString> coords = this->getCoordinatesIdList(name, path);
     std::list<StdString>::const_iterator it = coords.begin(), end = coords.end();
@@ -514 +517 @@
     {
       const StdString& coord = *it;
-      if (this->hasVariable(coord, path))
+      if (this->hasVariable(coord, path) && !this->isTemporal(coord, path))
       {
         std::map<StdString, StdSize> dimvar = this->getDimensions(&coord, path);
-        if (dimvar.size() != 2) return false;
-      }
-      else return false;
-    }
-    return true;
+        if (2 == dimvar.size()) ++nbLonLat;
+      }
+    }
+    if (2 != nbLonLat) isCurVi = false;
+
+    return isCurVi;
   }
 
@@ -538 +542 @@
     {
       const StdString& coord = *it;
-      if (this->hasVariable(coord, path))
+      if (this->hasVariable(coord, path) && !this->isTemporal(coord, path))
       {
         std::map<StdString, StdSize> dimvar = this->getDimensions(&coord, path);
@@ -581 +585 @@
       else
       {
-        if (coord.compare(this->getUnlimitedDimensionName()) == 0)
+        StdString unlimitedDimName = this->getUnlimitedDimensionName();
+        if (coord.compare(0, unlimitedDimName.size(), unlimitedDimName) == 0)
           continue;
         i++;
@@ -612 +617 @@
         else
         {
-          if (coord.compare(this->getUnlimitedDimensionName()) == 0)
+          StdString unlimitedDimName = this->getUnlimitedDimensionName();
+          if (coord.compare(0, unlimitedDimName.size(), unlimitedDimName) == 0)
             continue;
           return false;
@@ -712 +718 @@
                           const CVarPath* const path, StdSize record);
 
-  template <>
-  void CINetCDF4::getData(CArray<double, 1>& data, const StdString& var,
-                          bool collective, StdSize record,
-                          const std::vector<StdSize>* start /*= NULL*/,
-                          const std::vector<StdSize>* count /*= NULL*/)
-  {
-    int varid = this->getVariable(var);
-
-    if (this->mpi && collective)
-      CNetCdfInterface::varParAccess(ncidp, varid, NC_COLLECTIVE);
-    else if (this->mpi && !collective)
-      CNetCdfInterface::varParAccess(ncidp, varid, NC_INDEPENDENT);
-
-    std::vector<StdSize> sstart, scount;
-    StdSize array_size = 1;
-    this->getDataInfo(var, NULL, record, sstart, scount, array_size, start, count);
-
-    if (data.numElements() != array_size)
-    {
-      ERROR("CONetCDF4::getData(...)",
-            << "[ Array size = " << data.numElements()
-            << ", Data size = "  << array_size
-            << " ] Invalid array size");
-    }
-
-    CNetCdfInterface::getVaraType(ncidp, varid, &sstart[0], &scount[0], data.dataFirst());
-  }
-
   //---------------------------------------------------------------
 
@@ -745 +723 @@
                                        const CVarPath* const path)
   {
+    std::list<StdString>::const_iterator itbList, itList, iteList;
     std::list<StdString> clist = this->getCoordinatesIdList(varname, path);
-    if (this->hasCoordinates(varname, path))
-      return *clist.begin();
-    else
-      return *clist.rbegin();
+    itbList = clist.begin(); iteList = clist.end();
+    for (itList = itbList; itList != iteList; ++itList)
+    {
+      if (this->hasAttribute(CCFKeywords::XIOS_CF_units, &(*itList), path))
+      {
+        StdString unit = this->getAttributeValue(CCFKeywords::XIOS_CF_units, &(*itList), path);
+        if (CCFConvention::XIOS_CF_Longitude_units.end() != CCFConvention::XIOS_CF_Longitude_units.find(unit))
+          return *itList;
+      }
+    }
   }
 
@@ -755 +740 @@
                                        const CVarPath* const path)
   {
+    std::list<StdString>::const_iterator itbList, itList, iteList;
     std::list<StdString> clist = this->getCoordinatesIdList(varname, path);
-    if (this->hasCoordinates(varname, path))
-      return *(++clist.begin());
-    else
-      return *(++clist.rbegin());
+    itbList = clist.begin(); iteList = clist.end();
+    for (itList = itbList; itList != iteList; ++itList)
+    {
+      if (this->hasAttribute(CCFKeywords::XIOS_CF_units, &(*itList), path))
+      {
+        StdString unit = this->getAttributeValue(CCFKeywords::XIOS_CF_units, &(*itList), path);
+        if (CCFConvention::XIOS_CF_Latitude_units.end() != CCFConvention::XIOS_CF_Latitude_units.find(unit))
+          return *itList;
+      }
+    }
   }
 

XIOS/trunk/src/io/inetcdf4.hpp

@@ -98 +98 @@
                    StdSize record = UNLIMITED_DIM);
 
-      template <class T>
-      void getData(CArray<T, 1>& data, const StdString& var,
+      template <class T, int Ndim>
+      void getData(CArray<T, Ndim>& data, const StdString& var,
                    bool collective, StdSize record,
                    const std::vector<StdSize>* start = NULL,

XIOS/trunk/src/io/nc4_data_input.cpp

@@ -202 +202 @@
     \param [in] dimSizeMap Dimensions and and their corresponding names and size read from file
     \param [in] emelentPosition position of domain in grid
+    \param [in] fieldId id (or name) associated with the grid
   */
   void CNc4DataInput::readDomainAttributeValueFromFile(CDomain* domain, std::map<StdString, StdSize>& dimSizeMap,
@@ -227 +228 @@
     else if (this->isCurvilinear(fieldId))
     {
-
+      int ni = domain->ni;
+      int nj = domain->nj;
+      std::vector<StdSize> nBeginLatLon(2), nSizeLatLon(2);
+      nBeginLatLon[0] = domain->jbegin.getValue(); nBeginLatLon[1] = domain->ibegin.getValue();
+      nSizeLatLon[0]  = nj; nSizeLatLon[1] = ni;
+
+      StdString latName = this->getLatCoordName(fieldId);
+      domain->latvalue_2d.resize(ni,nj);
+      readFieldVariableValue(domain->latvalue_2d, latName, nBeginLatLon, nSizeLatLon);
+      StdString lonName = this->getLonCoordName(fieldId);
+      domain->lonvalue_2d.resize(ni,nj);
+      readFieldVariableValue(domain->lonvalue_2d, lonName, nBeginLatLon, nSizeLatLon);
+
+      StdString boundsLatName = this->getBoundsId(latName);
+      if (0 == boundsLatName.compare(""))
+         ERROR("CNc4DataInput::readDomainAttributeValueFromFile(...)",
+              << "Field '" << fieldId << std::endl
+              << "Trying to read attributes from curvilinear grid."
+              << "Latitude variable " << latName << " does not have bounds.");
+      StdString boundsLonName = this->getBoundsId(lonName);
+      if (0 == boundsLonName.compare(""))
+         ERROR("CNc4DataInput::readDomainAttributeValueFromFile(...)",
+              << "Field '" << fieldId << std::endl
+              << "Trying to read attributes from curvilinear grid."
+              << "Longitude variable " << lonName << " does not have bounds.");
+
+      domain->nvertex.setValue(4);
+      std::vector<StdSize> nBeginBndsLatLon(3), nSizeBndsLatLon(3);
+      nBeginBndsLatLon[0] = domain->jbegin.getValue(); nSizeBndsLatLon[0] = nj;
+      nBeginBndsLatLon[1] = domain->ibegin.getValue(); nSizeBndsLatLon[1] = ni;
+      nBeginBndsLatLon[2] = 0; nSizeBndsLatLon[2] = 4;
+
+      domain->bounds_lat_2d.resize(4,ni,nj);
+      readFieldVariableValue(domain->bounds_lat_2d, boundsLatName, nBeginBndsLatLon, nSizeBndsLatLon);
+      domain->bounds_lon_2d.resize(4,ni,nj);
+      readFieldVariableValue(domain->bounds_lon_2d, boundsLonName, nBeginBndsLatLon, nSizeBndsLatLon);
     }
     else if (this->isUnstructured(fieldId))
@@ -236 +272 @@
 
   /*!
-    Read attributes of a domain from a file
+    Read attribute value of a domain from a file
     \param [in] domain domain whose attributes are read from the file
     \param [in] dimSizeMap Dimensions and and their corresponding names and size read from file
     \param [in] emelentPosition position of domain in grid
+    \param [in] fieldId id (or name) associated with the grid
   */
   void CNc4DataInput::readDomainAttributesFromFile(CDomain* domain, std::map<StdString, StdSize>& dimSizeMap,
@@ -251 +288 @@
     itMapNi = itMapNj; ++itMapNi;
 
-    if (this->isRectilinear(fieldId))
+    if (this->isRectilinear(fieldId) || this->isCurvilinear(fieldId))
     {
       domain->nj_glo.setValue(itMapNj->second);
       domain->ni_glo.setValue((itMapNi)->second);
-    }
-    else if (this->isCurvilinear(fieldId))
-    {
-
     }
     else if (this->isUnstructured(fieldId))
@@ -271 +304 @@
     \param [in] dimSizeMap Dimensions and and their corresponding names and size read from file
     \param [in] emelentPosition position of axis in grid
+    \param [in] fieldId id (or name) associated with the grid
   */
   void CNc4DataInput::readAxisAttributesFromFile(CAxis* axis, std::map<StdString, StdSize>& dimSizeMap,
@@ -282 +316 @@
 
   /*!
-    Read attributes of an axis from a file
+    Read attribute value of an axis from a file
     \param [in] axis axis whose attributes are read from the file
     \param [in] dimSizeMap Dimensions and and their corresponding names and size read from file
     \param [in] emelentPosition position of axis in grid
+    \param [in] fieldId id (or name) associated with the grid
   */
   void CNc4DataInput::readAxisAttributeValueFromFile(CAxis* axis, std::map<StdString, StdSize>& dimSizeMap,
@@ -306 +341 @@
   }
 
-  void CNc4DataInput::readFieldVariableValue(CArray<double,1>& var, const StdString& varId,
-                                             const std::vector<StdSize>& nBegin,
-                                             const std::vector<StdSize>& nSize,
-                                             bool forceIndependent)
-  {
-    if (SuperClass::type==MULTI_FILE || !isCollective) return;
-
-    bool openCollective = isCollective;
-    if (forceIndependent) openCollective = !isCollective;
-    switch (SuperClass::type)
-    {
-      case MULTI_FILE:
-        SuperClassWriter::getData(var, varId, openCollective, 0);
-        break;
-      case ONE_FILE:
-      {
-        SuperClassWriter::getData(var, varId, openCollective, 0, &nBegin, &nSize);
-        break;
-      }
-    }
-  }
-
   void CNc4DataInput::closeFile_(void)
   {

XIOS/trunk/src/io/nc4_data_input.hpp

@@ -50 +50 @@
                                         int elementPosition, const StdString& fieldId);
 
-    void readFieldVariableValue(CArray<double,1>& var, const StdString& varId,
+    template <typename T, int Ndim>
+    void readFieldVariableValue(CArray<T, Ndim>& var, const StdString& varId,
                                 const std::vector<StdSize>& nBegin,
                                 const std::vector<StdSize>& nSize,
@@ -65 +66 @@
     bool isCollective;
   }; // class CNc4DataInput
+
+template <typename T, int Ndim>
+void CNc4DataInput::readFieldVariableValue(CArray<T,Ndim>& var, const StdString& varId,
+                                           const std::vector<StdSize>& nBegin,
+                                           const std::vector<StdSize>& nSize,
+                                           bool forceIndependent)
+{
+  if (SuperClass::type==MULTI_FILE || !isCollective) return;
+
+  bool openCollective = isCollective;
+  if (forceIndependent) openCollective = !isCollective;
+  switch (SuperClass::type)
+  {
+    case MULTI_FILE:
+      ERROR("CINetCDF4::readFieldVariableValue(...)",
+            << "Only support attributes reading with one_file mode");
+      break;
+    case ONE_FILE:
+    {
+      SuperClassWriter::getData(var, varId, openCollective, 0, &nBegin, &nSize);
+      break;
+    }
+  }
+}
 } // namespace xios
 

XIOS/trunk/src/io/netCdfInterface.cpp

@@ -603 +603 @@
 int CNetCdfInterface::inqAttName(int ncid, int varid, int attnum, StdString& name)
 {
-  char attName[NC_MAX_NAME + 1];
-  int status = nc_inq_attname(ncid, varid, attnum, attName);
+  std::vector<char> attName(NC_MAX_NAME + 1,' ');
+  int status = nc_inq_attname(ncid, varid, attnum, &attName[0]);
   if (NC_NOERR != status)
   {
@@ -616 +616 @@
     throw CNetCdfException(e);
   }
+
+  int nameSize = 0;
+  while ((nameSize < NC_MAX_NAME) && (' ' != attName[nameSize] )) ++nameSize;
+  name.resize(nameSize);
+//  for (int idx = 0; idx < nameSize; ++idx) name.at(idx) = attName[idx];
+  std::copy(&attName[0], &attName[nameSize-1], name.begin());
 
   return status;

XIOS/trunk/src/node/domain.cpp

@@ -228 +228 @@
    {
      if (this->isRedistributed_) return;
-     if (type_attr::rectilinear == type)
+     if ((type_attr::rectilinear == type)  || (type_attr::curvilinear == type))
      {
         this->isRedistributed_ = true;
@@ -328 +328 @@
 
         // Now fill other attributes
-        fillInRectilinearLonLat();
+        if (type_attr::rectilinear == type) fillInRectilinearLonLat();
      }
    }

XIOS/trunk/src/test/test_new_features.f90

@@ -27 +27 @@
 
   DOUBLE PRECISION,DIMENSION(ni_glo,nj_glo) :: lon_glo,lat_glo
+  DOUBLE PRECISION,DIMENSION(4,ni_glo,nj_glo) :: bnds_lon_glo, bnds_lat_glo
   DOUBLE PRECISION :: field_A_glo(ni_glo,nj_glo,llm), lval_ni_glo(ni_glo), lval_nj_glo(nj_glo)
   DOUBLE PRECISION,ALLOCATABLE :: lon(:,:),lat(:,:),field_A(:,:,:), field_All_Axis(:,:,:), lonvalue(:) , &
                                   field_Axis(:), lvaln(:), lval_ni(:), lval_nj(:), field_Two_Axis(:,:), lvalnInterp(:), &
-                                  lontransformed(:,:), lattransformed(:,:), lon_glotransformed(:,:), lat_glotransformed(:,:)
+                                  lontransformed(:,:), lattransformed(:,:), lon_glotransformed(:,:), lat_glotransformed(:,:), &
+                                  bnds_lon(:,:,:), bnds_lat(:,:,:)
   INTEGER :: ni,ibegin,iend,nj,jbegin,jend, nAxis, axisBegin, axisEnd
   INTEGER :: axisterpBegin, nAxisinterp, axisinterpEnd
   INTEGER :: niDomInterp,ibeginDomInterp,iendDomInterp,njDomInterp,jbeginDomInterp,jendDomInterp, niDomGlo, njDomGlo
-  INTEGER :: i,j,l,ts,n
+  INTEGER :: i,j,k,l,ts,n
 
 !!! MPI Initialization
@@ -53 +55 @@
       lon_glo(i,j)=(i-1)+(j-1)*ni_glo
       lat_glo(i,j)=1000+(i-1)+(j-1)*ni_glo
+      DO k = 1,4
+        bnds_lon_glo(1,i,j) = lon_glo(i,j)
+        bnds_lon_glo(2,i,j) = lon_glo(i,j)-10
+        bnds_lon_glo(3,i,j) = lon_glo(i,j)
+        bnds_lon_glo(4,i,j) = lon_glo(i,j)+10
+
+        bnds_lat_glo(1,i,j) = lat_glo(i,j) -10
+        bnds_lat_glo(2,i,j) = lat_glo(i,j)
+        bnds_lat_glo(3,i,j) = lat_glo(i,j) -10
+        bnds_lat_glo(4,i,j) = lat_glo(i,j)
+      ENDDO
       lval_ni_glo(i) = i-1
       DO l=1,llm
@@ -96 +109 @@
           lon(ni,nj),lat(ni,nj), lonvalue(ni*nj), &
           lvaln(nAxis), lval_ni(ni), lval_nj(nj), lvalnInterp(nAxisinterp), &
-          lontransformed(niDomInterp, njDomInterp), lattransformed(niDomInterp, njDomInterp))
+          lontransformed(niDomInterp, njDomInterp), lattransformed(niDomInterp, njDomInterp), &
+          bnds_lon(4,ni,nj), bnds_lat(4,ni,nj))
 
   ALLOCATE(mask(nj))
@@ -108 +122 @@
   lon(:,:)=lon_glo(ibegin+1:iend+1,jbegin+1:jend+1)
   lat(:,:)=lat_glo(ibegin+1:iend+1,jbegin+1:jend+1)
+  bnds_lon(:,:,:) = bnds_lon_glo(:,ibegin+1:iend+1,jbegin+1:jend+1)
+  bnds_lat(:,:,:) = bnds_lat_glo(:,ibegin+1:iend+1,jbegin+1:jend+1)
   lontransformed(:,:) = lon_glotransformed(ibeginDomInterp+1:iendDomInterp+1,jbeginDomInterp+1:jendDomInterp+1)
   lattransformed(:,:) = lat_glotransformed(ibeginDomInterp+1:iendDomInterp+1,jbeginDomInterp+1:jendDomInterp+1)
@@ -135 +151 @@
   CALL xios_set_domain_attr("domain_A",data_dim=2, data_ibegin=-1, data_ni=ni+2, data_jbegin=-2, data_nj=nj+4)
   CALL xios_set_domain_attr("domain_A",lonvalue_2D=lon,latvalue_2D=lat, type='curvilinear')
+  CALL xios_set_domain_attr("domain_A",bounds_lon_2d=bnds_lon,bounds_lat_2d=bnds_lat, nvertex=4, type='curvilinear')
 
   CALL xios_set_domain_attr("domain_A_transformed", ni_glo=niDomGlo, nj_glo=njDomGlo, &
@@ -201 +218 @@
   DEALLOCATE(field_A, field_Two_Axis, field_Axis, field_All_Axis, &
              lon, lat, lonvalue, &
-             lvaln, lval_ni, lval_nj, lvalnInterp)
+             lvaln, lval_ni, lval_nj, lvalnInterp, &
+             bnds_lon, bnds_lat)
 
   DEALLOCATE(mask)