![]() Saving this to paraview is as always: julia> Write_Paraview(GPS_Sanchez_grid, "GPSAlps_Sanchez_2017_grid") Also note that we do not attach units to the vector fields, but we do have them for the scalar fields: julia> GPS_Sanchez_grid = GeoData(lon,lat,height,(Velocity_mm_year=(Ve,Vn,Vz),V_north=Vn*mm/yr, V_east=Ve*mm/yr, V_vertical=Vz*mm/yr, Vmagnitude = Vmagnitude*mm/yr, Topography = height*km))įields: (:Velocity_mm_year, :V_north, :V_east, :V_vertical, :Vmagnitude, :Topography) As this transformation does not retain the east/north components of the velocity field, it is a good idea to save them as separate fields so we can color the vectors accordingly in Paraview. ![]() When saving to paraview, GMG internally does a vector transformation. As the velocity is a vector field, we need to save it as a data structure with 3 components. Saving and plotting in ParaviewĪt this stage, we have all we need. Julia> height = interpol.(lon,lat)/1e3 5. Julia> interpol = LinearInterpolation((Elevation.x, Elevation.y), Elevation.z') Next, we use the Interpolations.jl package to interpolate the topography: julia> using Interpolations We are using the GMT.jl to load the topographic data: julia> using GMT Yet, a better way is to load the topographic map of the area and interpolate the elevation to the velocity grid. We could ignore that and set the elevation to zero, which would allow saving the data directly to paraview. Yet, we don't have information yet about the elevation of the stations (as the provided data set did not give this). ![]() Interpolate topography on gridĪt this stage we have the 3D velocity components on a grid. Yet, given the small velocities in the Alps, it makes more sense to have them in units of mm/yr: julia> Vz = Vz*1000 Īnd the magnitude is: julia> Vmagnitude = sqrt.(Ve.^2 + Vn.^2 + Vz.^2) 4. Ind = intersect(findall(x->x=lon_Hz, lon), findall(x->x=lat_Hz, lat))Īt this stage, we have horizontal and vertical velocities in units of m/yr. Next we loop over all points in lon_Hz,lat_Hz and place them into the 2D matrixes: julia> for i in eachindex(lon_Hz) The strategy we take is to first define 2D matrixes with horizontal velocities with the same size as Vz which are initialized with NaN (not a number), which is treated specially by Paraview. So it appears that the horizontal velocities are given on the same regular grid as well, but not in the water. Let's plot the data as well: julia> Plots.scatter(lon_Hz,lat_Hz) Julia> lon_Hz, lat_Hz, Ve_Hz, Vn_Hz = data, data, data, data Julia> data = ParseColumns_CSV_File(data_file, 6) Next, we load the horizontal velocities from the file ALPS2017_DEF_HZ.GRD julia> data_file = CSV.File("ALPS2017_DEF_HZ.GRD",datarow=18,header=false,delim=' ') That is why we first initialize 3D matrixes for lon,lat,Vz: julia> lon, lat, Vz = zeros(41,31,1),zeros(41,31,1),zeros(41,31,1)Īnd we can reshape the vectors accordingly: julia> lon = reshape(lon_Vz,(41,31)) GMG requires 3D matrixes for the data (as we want to plot the results in paraview in 3D). ![]() We can determine the size of the grid with julia> unique(lon_Vz) So clearly, this is a fully regular grid. Let's have a look at the data, by plotting it: julia> using Plots Julia> lon_Vz, lat_Vz, Vz_vec = data, data, data 2. We can read the numerical data from the file with: julia> data = ParseColumns_CSV_File(data_file, 4) Julia> data_file = CSV.File("ALPS2017_DEF_VT.GRD",datarow=18,header=false,delim=' ') We can load that in julia as: julia> using CSV, GeophysicalModelGenerator So we have 4 columns with data values, and the data is separated by spaces. If we open it with a text editor, we see that the data starts at line 18, and has the following format: Column 1: Longitude Ĭolumn 3: Velocity in the height direction Ĭolumn 4: Uncertainty of the height component Let's have a look at the file ALPS2017_DEF_VT.GRD. We will use the package CSV.jl to load the comma-separated data. Next, we have a look at the data themselves. ALPS2017 DEFVT Vertical deformation model of the Alpine Region 2018/ALPS2017DEF_VT.GRD.ALPS2017 DEFHZ Surface deformation model of the Alpine Region 2018/ALPS2017DEF_HZ.GRD.Here, we will use the gridded data as an example (which interpolates the ), and will therefore download the following data sets: Some are the data on the actual stations and some are interpolated data on a grid. The data related to the paper can be downloaded from: There you will find links to several data sets. The example is based on a paper by Sanchez et al. The aim of this tutorial is to show you how you can download and plot GPS data, which are vector data.
0 Comments
Leave a Reply. |