BugFix

Processing/Create_HMSL-P06_vp.cpp

@@ -12,7 +12,7 @@ using namespace std;
 
 int main(){
 
-    auto tomo=Tomography::Model("HMSL-P06_dvp.nc");
+    auto tomo=Tomography("HMSL-P06_dvp.nc");
 
     // write nc file.
     int retval,ncid;
@@ -105,7 +105,7 @@ int main(){
     for (size_t i=0;i<tomo.GetDepths().size();++i) {
         for (size_t j=0;j<tomo.GetLatitudes().size();++j) {
             for (size_t k=0;k<tomo.GetLongitudes().size();++k) {
-                data[index]=Dvp_ak135(tomo.GetDepths()[i])*(1-tomo.GetValues()[index]/100.0);
+                data[index]=Dvp_ak135(tomo.GetDepths()[i])*(1+tomo.GetValues()[index]/100.0);
                 ++index;
             }
         }

Processing/Create_HMSL-S06_vs.cpp

@@ -12,7 +12,7 @@ using namespace std;
 
 int main(){
 
-    auto tomo=Tomography::Model("HMSL-S06_dvs.nc");
+    auto tomo=Tomography("HMSL-S06_dvs.nc");
 
     // write nc file.
     int retval,ncid;
@@ -105,7 +105,7 @@ int main(){
     for (size_t i=0;i<tomo.GetDepths().size();++i) {
         for (size_t j=0;j<tomo.GetLatitudes().size();++j) {
             for (size_t k=0;k<tomo.GetLongitudes().size();++k) {
-                data[index]=Dvs_ak135(tomo.GetDepths()[i])*(1-tomo.GetValues()[index]/100.0);
+                data[index]=Dvs_ak135(tomo.GetDepths()[i])*(1+tomo.GetValues()[index]/100.0);
                 ++index;
             }
         }

Processing/Create_SGLOBE-rani_vsv_vsh.cpp

@@ -77,18 +77,14 @@ int main(){
     retval=nc_get_var_float(ncid,lon_varid,lon);
     if (retval!=0) throw std::runtime_error(nc_strerror(retval));
 
-cout << depth_len << " " << lat_len << " " << lon_len << endl;
-
     // get dvs.
 
     size_t dvs_len=depth_len*lat_len*lon_len;
 
-cout << "Here" << endl;
     int dvs_varid;
     retval=nc_inq_varid(ncid,"dvs",&dvs_varid);
     if (retval!=0) throw std::runtime_error(nc_strerror(retval));
 
-cout << "There" << endl;
     float *dvs=new float [dvs_len];
     float *vs=new float [dvs_len];
     retval=nc_get_var_float(ncid,dvs_varid,dvs);

Processing/Create_TX2000_dvs.cpp

@@ -0,0 +1,187 @@
+#include<iostream>
+#include<fstream>
+#include<vector>
+#include<cmath>
+#include<set>
+#include<algorithm>
+
+#include<PREM.hpp>
+
+extern "C" {
+// Current version: netcdf-4.6.1
+#include<netcdf.h>
+}
+
+// the downloaded data is named "TX2000_percent.nc"
+// TX2000_percent.nc has longitude range -3 ~ 357.
+// Will adjust it to 1 ~ 357 (lose -3 values).
+
+using namespace std;
+
+int main(){
+
+    string infile="TX2000_percent.nc";
+    string outfile="TX2000_dvs.nc";
+
+    // read nc file.
+    int retval,ncid;
+    retval=nc_open(infile.c_str(),NC_NOWRITE,&ncid);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+    // get depth length and content.
+    int depth_id;
+    size_t depth_len;
+    retval=nc_inq_dimid(ncid,"depth",&depth_id);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+    retval=nc_inq_dimlen(ncid,depth_id,&depth_len);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+    int depth_varid;
+    retval=nc_inq_varid(ncid,"depth",&depth_varid);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+    float *depth = new float [depth_len];
+    retval=nc_get_var_float(ncid,depth_varid,depth);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+
+    // get latitude length and content.
+    int lat_id;
+    size_t lat_len;
+    retval=nc_inq_dimid(ncid,"latitude",&lat_id);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+    retval=nc_inq_dimlen(ncid,lat_id,&lat_len);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+    int lat_varid;
+    retval=nc_inq_varid(ncid,"latitude",&lat_varid);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+    float *lat=new float [lat_len];
+    retval=nc_get_var_float(ncid,lat_varid,lat);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+cout << lat[0] << endl;
+
+
+    // get longitude length and content.
+    int lon_id;
+    size_t lon_len;
+    retval=nc_inq_dimid(ncid,"longitude",&lon_id);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+    retval=nc_inq_dimlen(ncid,lon_id,&lon_len);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+    int lon_varid;
+    retval=nc_inq_varid(ncid,"longitude",&lon_varid);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+    float *lon_raw=new float [lon_len];
+    retval=nc_get_var_float(ncid,lon_varid,lon_raw);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+    // get dvs.
+
+    size_t dvs_raw_len=depth_len*lat_len*lon_len;
+
+    int dvs_varid;
+    retval=nc_inq_varid(ncid,"dvs",&dvs_varid);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+    float *dvs_raw=new float [dvs_raw_len];
+    retval=nc_get_var_float(ncid,dvs_varid,dvs_raw);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+    // close infile.
+    retval=nc_close(ncid);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+
+    // write nc file.
+    retval=nc_create(outfile.c_str(),NC_CLOBBER,&ncid);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+
+    // define dimension id.
+    int depth_dimid,lat_dimid,lon_dimid;
+    retval=nc_def_dim(ncid,"depth",depth_len,&depth_dimid);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+    retval=nc_def_dim(ncid,"latitude",lat_len,&lat_dimid);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+    retval=nc_def_dim(ncid,"longitude",lon_len-1,&lon_dimid);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+    // v is not a dimension, just a variable,
+    // we will generate a variable id from the dimension ids from dimensions.
+    int dimids[3];
+    dimids[0]=depth_dimid;
+    dimids[1]=lat_dimid;
+    dimids[2]=lon_dimid;
+
+
+    // define variable id.
+    retval=nc_def_var(ncid,"depth",NC_FLOAT,1,&depth_dimid,&depth_varid);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+    retval=nc_def_var(ncid,"latitude",NC_FLOAT,1,&lat_dimid,&lat_varid);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+    retval=nc_def_var(ncid,"longitude",NC_FLOAT,1,&lon_dimid,&lon_varid);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+    int v_varid;
+    retval=nc_def_var(ncid,"v",NC_FLOAT,3,dimids,&v_varid);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+    // define units.
+    string unit_depth="km_downward",unit_lat="degrees_south",unit_lon="degrees_east";
+    retval=nc_put_att_text(ncid,depth_varid,"units",unit_depth.size(),unit_depth.c_str());
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+    retval=nc_put_att_text(ncid,lat_varid,"units",unit_lat.size(),unit_lat.c_str());
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+    retval=nc_put_att_text(ncid,lon_varid,"units",unit_lon.size(),unit_lon.c_str());
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+    string unit_v="Vs perturbation (\% deviation from layer mean)";
+    retval=nc_put_att_text(ncid,v_varid,"units",unit_v.size(),unit_v.c_str());
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+    // End of variable definition.
+    retval=nc_enddef(ncid);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+    // load data.
+
+    /// depth.
+    retval=nc_put_var_float(ncid,depth_varid,depth);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+    /// lat.
+cout << lat[0] << endl;
+    retval=nc_put_var_float(ncid,lat_varid,lat);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+    /// lon.
+
+    float *lon=new float [lon_len-1];
+    for (size_t i=1;i<lon_len;++i) lon[i-1]=lon_raw[i];
+    retval=nc_put_var_float(ncid,lon_varid,lon);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+    /// v.
+    size_t dvs_len=depth_len*lat_len*(lon_len-1);
+    float *dvs=new float [dvs_len];
+    for (size_t i=0;i<dvs_raw_len;++i) {
+        size_t depth_index=i/(lat_len*lon_len);
+        size_t lat_index=(i%(lat_len*lon_len))/lon_len;
+        size_t lon_index=(i%(lat_len*lon_len))%lon_len;
+        if (lon_index==0) continue;
+        --lon_index;
+        dvs[depth_index*(lat_len*(lon_len-1))+lat_index*(lon_len-1)+lon_index]=dvs_raw[i];
+    }
+    retval=nc_put_var_float(ncid,v_varid,dvs);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+    // close the nc file.
+    retval=nc_close(ncid);
+    if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+    return 0;
+}

Processing/Create_TX2011_vs.cpp

@@ -14,7 +14,7 @@ double TX2011_ref(const double &d);
 
 int main(){
 
-    auto tomo=Tomography::Model("TX2011_dvs.nc");
+    auto tomo=Tomography("TX2011_dvs.nc");
 
     // write nc file.
     int retval,ncid;
@@ -107,7 +107,7 @@ int main(){
     for (size_t i=0;i<tomo.GetDepths().size();++i) {
         for (size_t j=0;j<tomo.GetLatitudes().size();++j) {
             for (size_t k=0;k<tomo.GetLongitudes().size();++k) {
-                data[index]=TX2011_ref(tomo.GetDepths()[i])*(1-tomo.GetValues()[index]/100.0);
+                data[index]=TX2011_ref(tomo.GetDepths()[i])*(1+tomo.GetValues()[index]/100.0);
                 ++index;
             }
         }

src/Tomography.hpp

@@ -28,17 +28,19 @@ extern "C" {
 class Tomography {
 
     std::vector<double> depth,lat,lon,v; // depth in km. lon between [-180,180].
-    bool Use360;
+    bool range_is_360;
 
 public:
     Tomography()=default;
+
     Tomography(const std::vector<double> &d, const std::vector<double> &lo,
           const std::vector<double> &la, const std::vector<double> &data) {
         if (data.size()!=d.size()*lo.size()*la.size())
             throw std::runtime_error("In "+std::string(__func__)+", input sizes don't match.");
-        depth=d;lat=la;lon=lo;v=data;Use360=(!lo.empty() && lo[0]>=0);
-
+        depth=d;lat=la;lon=lo;v=data;range_is_360=(!lon.empty() && lon[0]>=0);
+        if (!la.empty() && la[0]>0) ReverseLatitude();
     }
+
     Tomography(const std::string &nc_filename){
 
         // tmp data.
@@ -89,7 +91,6 @@ class Tomography {
         for (size_t i=0;i<lat_len;++i) lat[i]=data[i];
         delete [] data;
 
-
         // get longitude length and data.
         int lon_id;
         size_t lon_len;
@@ -106,10 +107,10 @@ class Tomography {
         data=new float [lon_len];
         retval=nc_get_var_float(ncid,lon_varid,data);
         if (retval!=0) throw std::runtime_error(nc_strerror(retval));
-        Use360=(data[0]>=0);
         for (size_t i=0;i<lon_len;++i) lon[i]=data[i];
         delete [] data;
 
+        range_is_360=(!lon.empty() && lon[0]>=0);
 
         // get velocity length and data.
         size_t v_len=depth_len*lat_len*lon_len;
@@ -129,13 +130,17 @@ class Tomography {
         // close the nc file.
         retval=nc_close(ncid);
         if (retval!=0) throw std::runtime_error(nc_strerror(retval));
+
+        // check latitude direction.
+        if (!lat.empty() && lat[0]>0) ReverseLatitude();
     }
 
     double GetValueAt(const double &d, double lo, const double &la) const;
     const std::vector<double> &GetDepths() const {return depth;}
     const std::vector<double> &GetLatitudes() const {return lat;}
     const std::vector<double> &GetLongitudes() const {return lon;}
     const std::vector<double> &GetValues() const {return v;}
+    void ReverseLatitude();
 
 };
 
@@ -146,26 +151,29 @@ double Tomography::GetValueAt(const double &d, double lo, const double &la) cons
 
     if (d<depth[0] || d>depth.back() || la<lat[0] || la>lat.back()) return 0.0/0.0;
 
-    lo=(Use360?Lon2360(lo):Lon2180(lo));
+    size_t lat_len=GetLatitudes().size(),lon_len=GetLongitudes().size();
+
+    lo=(range_is_360?Lon2360(lo):Lon2180(lo));
     size_t index_depth=std::distance(depth.begin(),
                                      std::lower_bound(depth.begin(),depth.end(),d));
     size_t index_lat=std::distance(lat.begin(),std::lower_bound(lat.begin(),lat.end(),la));
     size_t index_lon=std::distance(lon.begin(),std::lower_bound(lon.begin(),lon.end(),lo));
 
     size_t index_prev_depth=(index_depth==0?0:index_depth-1);
     size_t index_prev_lat=(index_lat==0?0:index_lat-1);
-    size_t index_prev_lon=(index_lon==0?lon.size()-1:index_lon-1);
-    if (index_lon==lon.size()) index_lon=0;
+    size_t index_prev_lon=(index_lon==0?lon_len-1:index_lon-1);
+    if (index_lon==lon_len) index_lon=0;
+
 
     std::vector<double> p(8),pp(3);
-    p[0]=v[index_prev_depth*lat.size()*lon.size()+index_prev_lat*lon.size()+index_prev_lon];
-    p[1]=v[index_prev_depth*lat.size()*lon.size()+index_prev_lat*lon.size()+index_lon];
-    p[2]=v[index_prev_depth*lat.size()*lon.size()+index_lat*lon.size()+index_prev_lon];
-    p[3]=v[index_prev_depth*lat.size()*lon.size()+index_lat*lon.size()+index_lon];
-    p[4]=v[index_depth*lat.size()*lon.size()+index_prev_lat*lon.size()+index_prev_lon];
-    p[5]=v[index_depth*lat.size()*lon.size()+index_prev_lat*lon.size()+index_lon];
-    p[6]=v[index_depth*lat.size()*lon.size()+index_lat*lon.size()+index_prev_lon];
-    p[7]=v[index_depth*lat.size()*lon.size()+index_lat*lon.size()+index_lon];
+    p[0]=v[index_prev_depth*lat_len*lon_len+index_prev_lat*lon_len+index_prev_lon];
+    p[1]=v[index_prev_depth*lat_len*lon_len+index_prev_lat*lon_len+index_lon];
+    p[2]=v[index_prev_depth*lat_len*lon_len+index_lat*lon_len+index_prev_lon];
+    p[3]=v[index_prev_depth*lat_len*lon_len+index_lat*lon_len+index_lon];
+    p[4]=v[index_depth*lat_len*lon_len+index_prev_lat*lon_len+index_prev_lon];
+    p[5]=v[index_depth*lat_len*lon_len+index_prev_lat*lon_len+index_lon];
+    p[6]=v[index_depth*lat_len*lon_len+index_lat*lon_len+index_prev_lon];
+    p[7]=v[index_depth*lat_len*lon_len+index_lat*lon_len+index_lon];
     if (index_prev_lon>index_lon)
         pp[0]=((lo<lon[index_prev_lon]?lo+360:lo)-lon[index_prev_lon])
               /(lon[index_lon]+360-lon[index_prev_lon]);
@@ -176,4 +184,18 @@ double Tomography::GetValueAt(const double &d, double lo, const double &la) cons
     return LinearInterp(p,pp);
 }
 
+void Tomography::ReverseLatitude(){
+    auto old_data=GetValues();
+    size_t lat_len=GetLatitudes().size(),lon_len=GetLongitudes().size();
+
+    for (size_t i=0;i<GetValues().size();++i){
+        size_t depth_index=i/(lat_len*lon_len);
+        size_t lat_index=lat_len-1-(i%(lat_len*lon_len))/lon_len;
+        size_t lon_index=(i%(lat_len*lon_len))%lon_len;
+        v[i]=old_data[depth_index*(lat_len*lon_len)+lat_index*lon_len+lon_index];
+    }
+    std::reverse(lat.begin(),lat.end());
+    return;
+}
+
 #endif