GPLIB++
ReadWriteImpedances.cpp
Go to the documentation of this file.
1 //============================================================================
2 // Name : ReadWriteImpedance.cpp
3 // Author : Jul 13, 2009
4 // Version :
5 // Copyright : 2009, mmoorkamp
6 //============================================================================
7 
8 #include <netcdfcpp.h>
9 #include <boost/bind.hpp>
10 #include <boost/function.hpp>
11 #include "../Global/Util.h"
12 #include "ReadWriteImpedances.h"
13 
14 namespace gplib
15  {
16  //! The name used for the x-coordinates of the measurements in netcdf files
17  static const std::string MeasPosXName = "MeasPosX";
18  //! The name used for the y-coordinates of the measurements in netcdf files
19  static const std::string MeasPosYName = "MeasPosY";
20  //! The name used for the z-coordinates of the measurements in netcdf files
21  static const std::string MeasPosZName = "MeasPosZ";
22  //! The name used for the index of the measurements in netcdf files
23  static const std::string StationNumberName = "StationNumber";
24 
25  //! Read a vector from a netcdf file
26  template<class VectorType>
27  void ReadVec(NcFile &NetCDFFile, const std::string &DataName,
28  const std::string &DimName, VectorType &Position)
29  {
30  //create a netcdf dimension for the Station number
31  NcDim *Dim = NetCDFFile.get_dim(DimName.c_str());
32  //determine the size of that dimension
33  const size_t nvalues = Dim->size();
34 
35  //allocate memory in the class variable
36  Position.resize(nvalues);
37  // create netcdf variable with the same name as the dimension
38  NcVar *SizeVar = NetCDFFile.get_var(DataName.c_str());
39  //read coordinate values from netcdf file
40  SizeVar->get(&Position[0], nvalues);
41  }
42 
43  //! Write a vectorial quantity to a netcdf file
44  template<class VectorType>
45  void WriteVec(NcFile &NetCDFFile, const std::string &MeasPosName,
46  const VectorType &Position, NcDim *Dimension, const std::string unit)
47  {
48  const size_t nmeas = Position.size();
49  NcVar *PosVar = NetCDFFile.add_var(MeasPosName.c_str(), ncDouble,
50  Dimension);
51  PosVar->add_att("units", unit.c_str());
52  PosVar->put(&Position[0], nmeas);
53  }
54 
55  const std::string FreqDimName = "Frequency";
56  //write one compoment of the impedance tensor to a netcdf file
57  //this is an internal helper function
58  void WriteImpedanceComp(NcFile &NetCDFFile, NcDim *StatNumDim,
59  NcDim *FreqDim, const gplib::rvec &Impedances,
60  const std::string &CompName, const size_t compindex)
61  {
62  NcVar *CompVar = NetCDFFile.add_var(CompName.c_str(), ncDouble,
63  FreqDim, StatNumDim);
64  gplib::rvec Component(FreqDim->size() * StatNumDim->size());
65  for (size_t i = 0; i < Component.size(); ++i)
66  {
67  Component(i) = Impedances(i * 8 + compindex);
68  }
69  CompVar->put(&Component[0], FreqDim->size(), StatNumDim->size());
70  }
71  //read one compoment of the impedance tensor from a netcdf file
72  //this is an internal helper function
73  void ReadImpedanceComp(NcFile &NetCDFFile, gplib::rvec &Impedances,
74  const std::string &CompName, const size_t compindex)
75  {
76  NcVar *SizeVar = NetCDFFile.get_var(CompName.c_str());
77  const size_t nvalues = Impedances.size() / 8;
78  gplib::rvec Temp(nvalues);
79  SizeVar->get(&Temp[0], SizeVar->edges()[0], SizeVar->edges()[1]);
80  for (size_t i = 0; i < nvalues; ++i)
81  {
82  Impedances(i * 8 + compindex) = Temp(i);
83  }
84 
85  }
86 
87  void WriteImpedancesToNetCDF(const std::string &filename,
88  const std::vector<double> &Frequencies,
89  const std::vector<double> &StatXCoord,
90  const std::vector<double> &StatYCoord,
91  const std::vector<double> &StatZCoord, const gplib::rvec &Impedances)
92  {
93  const size_t nstats = StatXCoord.size();
94  const size_t nfreqs = Frequencies.size();
95  assert(nstats == StatYCoord.size());
96  assert(nstats == StatYCoord.size());
97  assert(Impedances.size() == nstats * nfreqs * 8);
98  //create a netcdf file
99  NcFile DataFile(filename.c_str(), NcFile::Replace);
100  //each station gets an index number that we use as a dimension
101  //in the netcdf file
102  NcDim *StatNumDim = DataFile.add_dim(StationNumberName.c_str(), nstats);
103  std::vector<int> StationNumber;
104  std::generate_n(back_inserter(StationNumber), nstats, IntSequence(0));
105  NcVar *StatNumVar = DataFile.add_var(StationNumberName.c_str(), ncInt,
106  StatNumDim);
107  StatNumVar->put(&StationNumber[0], nstats);
108  //write out the measurement coordinates
109  WriteVec(DataFile, MeasPosXName, StatXCoord, StatNumDim, "m");
110  WriteVec(DataFile, MeasPosYName, StatYCoord, StatNumDim, "m");
111  WriteVec(DataFile, MeasPosZName, StatZCoord, StatNumDim, "m");
112  //write out the frequencies that we store
113  NcDim *FreqDim = DataFile.add_dim(FreqDimName.c_str(),
114  Frequencies.size());
115  NcVar *FreqVar = DataFile.add_var(FreqDimName.c_str(), ncDouble,
116  FreqDim);
117  FreqVar->put(&Frequencies[0], nfreqs);
118  //and now we can write all the impedance components
119  WriteImpedanceComp(DataFile, StatNumDim, FreqDim, Impedances, "Zxx_re",
120  0);
121  WriteImpedanceComp(DataFile, StatNumDim, FreqDim, Impedances, "Zxx_im",
122  1);
123  WriteImpedanceComp(DataFile, StatNumDim, FreqDim, Impedances, "Zxy_re",
124  2);
125  WriteImpedanceComp(DataFile, StatNumDim, FreqDim, Impedances, "Zxy_im",
126  3);
127  WriteImpedanceComp(DataFile, StatNumDim, FreqDim, Impedances, "Zyx_re",
128  4);
129  WriteImpedanceComp(DataFile, StatNumDim, FreqDim, Impedances, "Zyx_im",
130  5);
131  WriteImpedanceComp(DataFile, StatNumDim, FreqDim, Impedances, "Zyy_re",
132  6);
133  WriteImpedanceComp(DataFile, StatNumDim, FreqDim, Impedances, "Zyy_im",
134  7);
135  }
136 
137  void ReadImpedancesFromNetCDF(const std::string &filename, std::vector<
138  double> &Frequencies, std::vector<double> &StatXCoord, std::vector<
139  double> &StatYCoord, std::vector<double> &StatZCoord,
140  gplib::rvec &Impedances)
141  {
142  NcFile DataFile(filename.c_str(), NcFile::ReadOnly);
143  ReadVec(DataFile, MeasPosXName, StationNumberName, StatXCoord);
144  ReadVec(DataFile, MeasPosYName, StationNumberName, StatYCoord);
145  ReadVec(DataFile, MeasPosZName, StationNumberName, StatZCoord);
146  ReadVec(DataFile, FreqDimName, FreqDimName, Frequencies);
147  Impedances.resize(Frequencies.size() * StatXCoord.size() * 8);
148  ReadImpedanceComp(DataFile, Impedances, "Zxx_re", 0);
149  ReadImpedanceComp(DataFile, Impedances, "Zxx_im", 1);
150  ReadImpedanceComp(DataFile, Impedances, "Zxy_re", 2);
151  ReadImpedanceComp(DataFile, Impedances, "Zxy_im", 3);
152  ReadImpedanceComp(DataFile, Impedances, "Zyx_re", 4);
153  ReadImpedanceComp(DataFile, Impedances, "Zyx_im", 5);
154  ReadImpedanceComp(DataFile, Impedances, "Zyy_re", 6);
155  ReadImpedanceComp(DataFile, Impedances, "Zyy_im", 7);
156 
157  }
158  }
void ReadVec(NcFile &NetCDFFile, const std::string &DataName, const std::string &DimName, VectorType &Position)
Read a vector from a netcdf file.
void WriteVec(NcFile &NetCDFFile, const std::string &MeasPosName, const VectorType &Position, NcDim *Dimension, const std::string unit)
Write a vectorial quantity to a netcdf file.
void WriteImpedanceComp(NcFile &NetCDFFile, NcDim *StatNumDim, NcDim *FreqDim, const gplib::rvec &Impedances, const std::string &CompName, const size_t compindex)
void ReadImpedanceComp(NcFile &NetCDFFile, gplib::rvec &Impedances, const std::string &CompName, const size_t compindex)
const std::string FreqDimName
void ReadImpedancesFromNetCDF(const std::string &filename, std::vector< double > &Frequencies, std::vector< double > &StatXCoord, std::vector< double > &StatYCoord, std::vector< double > &StatZCoord, gplib::rvec &Impedances)
Read magnetotelluric impedances from a netcdf file.
void WriteImpedancesToNetCDF(const std::string &filename, const std::vector< double > &Frequencies, const std::vector< double > &StatXCoord, const std::vector< double > &StatYCoord, const std::vector< double > &StatZCoord, const gplib::rvec &Impedances)
Write magnetotelluric impedances to a netcdf file.