Hi all,
I am hoping that I could get some feedback both for my code and for the plots they’ve generated. If this isn’t the correct forum for this type of post, please let me know where a better suited place would be.
For some context, these plots are showing residual GPS values. This means, for those unfamiliar, that certain loading forces have been removed from this model, such as non-tidal atmospheric and oceanic loading, seasonal loading (snow for example), hydrological loading (large bodies of (non oceanic) water that is controlled for seasonality), among others.
As for the code, I am really interested in learning how to streamline and create more efficient scripts. I’d consider myself to be a novice still, and would love to continue to learn how to be more efficient in my implementation.
For my plots, there is a little less wiggle-room for change, since my grad advisor requested the shown attributed, but I am happy to make minor changes as long as the data being presented remains easy to read.
Here is my code, and I’m happy to provide some of my data if anyone is interested:
using DelimitedFiles, Glob, GMT
# paths to various directories
imgpath = "/home/rob/Pictures/gmtpng/stack_plot/";
twosig = "/home/rob/Documents/gmt_jl/twosig_ten/";
values = "/home/rob/Documents/gmt_jl/value_match/";
lonlat = gmtread("lonlat.txt");
## Function to read each file in as a GMT dataset ##
function gmt_read(i)
Sepoch = gmtread(i[1])
Vepoch = gmtread(i[2])
return Sepoch, Vepoch
end
for i in zip((glob("*.txt", twosig)), (glob("*.txt", values)))
Sepoch, Vepoch = gmt_read(i)
C = GMT.makecpt(cmap=:vik,
range=(-4, 4, 0.1)
);
G1 = GMT.gridit(Sepoch, mask = 0.3, inc = 0.1);
G2 = GMT.gridit(Vepoch, mask = 0.3, inc = 0.1);
base = basename(i[1])
yr = string(base[1:8])
GMT.subplot(grid = "2x1", savefig = imgpath * "$yr" * "_stack.png",);
# gridding image for 2sigma values
GMT.grdimage(G1, region = [-130 -70 24 51],
proj = (name = :lambertConic, center = [-100 35], parallels = [33 45]),
cmap = C,
colorbar = (pos = (outside = true, anchor = :MR, offset = (0, .9)), xlabel = "mm"),
coast = (borders = ((type = 1, pen = ("thick", "black")), (type = 2, pen = ("thinner", "lightgrey"))),
area = 500,
shore = (:thinnest, :lightgrey)),
panel = (1,1),
title = "Vertical Residual greater than 2sigma | $yr",
);
# gridding image for values
GMT.grdimage(G2, region = [-130 -70 24 51],
proj = (name = :lambertConic, center = [-100 35], parallels = [33 45]),
cmap = C,
colorbar = (pos = (outside = true, anchor = :MR, offset = (0, .9)), xlabel = "mm"),
coast = (borders = ((type = 1, pen = ("thick", "black")), (type = 2, pen = ("thinner", "lightgrey"))),
area = 500,
shore = (:thinnest, :lightgrey)),
panel = (2,1),
title = "Vertical Residual | $yr",
);
GMT.scatter!(lonlat,
marker = :c,
markersize = 0.075,
mc = "black",
#figname = imgpath * "$yr" *" _stacked.png",
panel = (1,1),
# show = true
);
GMT.scatter!(lonlat,
marker = :c,
markersize = 0.075,
mc = "black",
#figname = imgpath * "$yr" *" _stacked.png",
panel = (2,1),
# show = true
);
subplot(:end)
end
Which produced plots that look like this (I’ll post 3 here, but this loops created 1,180 plots):
As of this post, the 2002.127 is the most recent plot, and I think it will take about 3-hours for Julia to create and save all my plots up to 2005.998.
Happy to clarify any confusion and answer any questions (to the best of my ability).
Thank you for any feedback and considerations!
- Rob
