Skip to content

Commit

Permalink
SP12RTS
Browse files Browse the repository at this point in the history
  • Loading branch information
@shuleyu
shuleyu committed Sep 15, 2018
1 parent dd2db21 commit d3bbcb2
Show file tree
Hide file tree
Showing 11 changed files with 903 additions and 3 deletions.
139 changes: 139 additions & 0 deletions Create_SP12RTS_dvp.cpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,139 @@
#include<iostream>
#include<fstream>
#include<vector>
#include<cmath>
#include<set>
#include<algorithm>

extern "C" {
// Current version: netcdf-4.6.1
#include<netcdf.h>
}

// the data generated by Create_SP12RTS_dvp.sh is named as "SP12RTS_dvp.txt";

using namespace std;

int main(){

vector<vector<float>> raw_data;
float lat,lon,depth,dvp;
set<float> s_lat,s_lon,s_depth;
string tmpstr;
ifstream fpin("SP12RTS_dvp.txt");

while (fpin >> depth >> lon >> lat >> dvp){
raw_data.push_back({depth,lat,lon,dvp});
s_lat.insert(lat);
s_lon.insert(lon);
s_depth.insert(depth);
}
fpin.close();

auto cmp=[](const vector<float> &v1, const vector<float> &v2){
if (v1[0]==v2[0]) {
if (v1[1]==v2[1]) return v1[2]<v2[2];
else return v1[1]<v2[1];
}
else return v1[0]<v2[0];
};
sort(raw_data.begin(),raw_data.end(),cmp);

// write nc file.
int retval,ncid;
string outfile="SP12RTS_dvp.nc";

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",s_depth.size(),&depth_dimid);
if (retval!=0) throw std::runtime_error(nc_strerror(retval));
retval=nc_def_dim(ncid,"latitude",s_lat.size(),&lat_dimid);
if (retval!=0) throw std::runtime_error(nc_strerror(retval));
retval=nc_def_dim(ncid,"longitude",s_lon.size(),&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.
int depth_varid,lat_varid,lon_varid;
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_north",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="dvp, (%), relative to PREM";
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.

float *data;
size_t i;

/// depth.
data=new float [s_depth.size()];
i=0;
for (auto &item:s_depth) data[i++]=item;
retval=nc_put_var_float(ncid,depth_varid,data);
if (retval!=0) throw std::runtime_error(nc_strerror(retval));
delete [] data;

/// lat.
data=new float [s_lat.size()];
i=0;
for (auto &item:s_lat) data[i++]=item;
retval=nc_put_var_float(ncid,lat_varid,data);
if (retval!=0) throw std::runtime_error(nc_strerror(retval));
delete [] data;

/// lon.
data=new float [s_lon.size()];
i=0;
for (auto &item:s_lon) data[i++]=item;
retval=nc_put_var_float(ncid,lon_varid,data);
if (retval!=0) throw std::runtime_error(nc_strerror(retval));
delete [] data;

/// v.
data=new float [raw_data.size()];
i=0;
for (auto &item:raw_data) data[i++]=item[3];
retval=nc_put_var_float(ncid,v_varid,data);
if (retval!=0) throw std::runtime_error(nc_strerror(retval));
delete [] data;

// close the nc file.
retval=nc_close(ncid);
if (retval!=0) throw std::runtime_error(nc_strerror(retval));

return 0;
}
159 changes: 159 additions & 0 deletions Create_SP12RTS_dvp.sh
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,159 @@
#!/bin/bash

# prompt> export TOMO_S20RTS=`pwd`
# prompt> export PATH=${PATH}:.
# SP12RTS, dvp (EP), the first argument to mkmap is 3.
# This will create SP12RTS_dvp.txt. Noitice it's not sorted.

rm -f SP12RTS_dvp.txt *xyz *ps
while read depth
do
./mkmap 3 1 ${depth}
Line=`tail -n 2 map.${depth}km.xyz`
NorthVal=`echo ${Line} | awk '{print $3}'`
SouthVal=`echo ${Line} | awk '{print $6}'`

for lon in `seq -180 179`
do
printf "%.2lf %.2lf %.2lf %.5e\n" ${depth} ${lon} -90 ${SouthVal} >> SP12RTS_dvp.txt
printf "%.2lf %.2lf %.2lf %.5e\n" ${depth} ${lon} 90 ${NorthVal} >> SP12RTS_dvp.txt
done
awk -v d=${depth} 'NR<=64440 {printf "%.2lf %.2lf %.2lf %.5e\n",d,$1,$2,$3}' map.${depth}km.xyz >> SP12RTS_dvp.txt

done << EOF
25
50
70
75
100
120
125
150
170
175
200
220
225
250
270
275
300
320
325
350
370
375
400
420
425
450
470
475
500
520
525
550
575
580
600
620
625
650
670
675
700
720
725
750
775
780
800
825
850
875
900
925
950
975
1000
1025
1050
1075
1100
1125
1150
1175
1200
1225
1250
1275
1300
1325
1350
1375
1400
1425
1450
1475
1500
1525
1550
1575
1600
1625
1650
1675
1700
1725
1750
1775
1800
1825
1850
1875
1900
1925
1950
1975
2000
2025
2050
2075
2100
2125
2150
2175
2200
2225
2250
2275
2300
2325
2350
2375
2400
2425
2450
2475
2500
2525
2550
2575
2600
2625
2650
2675
2700
2725
2750
2775
2800
2825
2850
2875
2891
EOF

rm -f *xyz *ps

exit 0
Loading

0 comments on commit d3bbcb2

Please sign in to comment.