use Literate

docs/Project.toml

@@ -2,7 +2,12 @@
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 GeoDatasets = "ddc7317b-88db-5cb5-a849-8449e5df04f9"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
-PythonPlot = "274fc56d-3b97-40fa-a1cd-1b4a50311bf9"
+PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 ROMS = "da6b2582-fac1-11ea-b7c4-a321ee60f94e"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[compat]
+Documenter = "1"
+Literate = "2"

docs/make.jl

@@ -1,13 +1,33 @@
 using Pkg
 Pkg.activate(@__DIR__)
 CI = get(ENV, "CI", nothing) == "true"
-using Documenter, ROMS
+using Documenter
+using ROMS
+import Literate
+
+Literate.markdown(
+    joinpath(@__DIR__, "..", "examples", "plots.jl"),
+    joinpath(@__DIR__, "src"),
+    execute = true,
+    documenter = true,
+    # We add the credit to Literate.jl the footer
+    credit = false,
+)
+
+Literate.notebook(
+    joinpath(@__DIR__, "..", "examples", "plots.jl"),
+    joinpath(@__DIR__, "src"),
+    execute = false,
+)
+
 
 makedocs(
     modules = [ROMS],
     sitename = "ROMS.jl",
     format = Documenter.HTML(
         prettyurls = CI,
+        footer = "Powered by [Documenter.jl](https://github.com/JuliaDocs/Documenter.jl), [Literate.jl](https://github.com/fredrikekre/Literate.jl) and the [Julia Programming Language](https://julialang.org/)"
+
     ),
     pages = [
         "Tutorial" => "index.md",
@@ -16,7 +36,5 @@ makedocs(
     ],
 )
 
-if CI
-    deploydocs(repo = "github.com/Alexander-Barth/ROMS.jl.git",
-               target = "build")
-end
+deploydocs(repo = "github.com/Alexander-Barth/ROMS.jl.git",
+           target = "build")

examples/plots.jl

@@ -0,0 +1,231 @@
+# # Plotting ROMS results and input files
+
+#md # The code here is also available as a [notebook](plots.ipynb).
+
+# The aim here is to visualize the model files with generic plotting and analsis packages rather than to use a model specific visualization tool which hides many details and might lack of flexibility.
+# The necessary files are already in the directory containing the model simulation and its
+# parent direction (`ROMS-implementation-test`). Downloading the files is only needed if you did not run the simulation.
+
+grid_fname = "LS2v.nc"
+
+if !isfile(grid_fname)
+    download("https://dox.ulg.ac.be/index.php/s/J9DXhUPXbyLADJa/download",grid_fname)
+end
+
+fname = "roms_his.nc"
+if !isfile(fname)
+    download("https://dox.ulg.ac.be/index.php/s/17UWsY7tRNMDf4w/download",fname)
+end
+
+# ## Bathymetry
+
+# In this example, the bathymetry defined in the grid file is visualized. Make sure that your current working directory
+# contains the file `LS2v.nc` (use e.g. `;cd ~/ROMS-implementation-test`)
+
+using ROMS, PyPlot, NCDatasets, GeoDatasets, Statistics
+
+ds_grid = NCDataset("LS2v.nc");
+lon = ds_grid["lon_rho"][:,:];
+lat = ds_grid["lat_rho"][:,:];
+h = nomissing(ds_grid["h"][:,:],NaN);
+mask_rho = ds_grid["mask_rho"][:,:];
+
+figure(figsize=(7,4))
+hmask = copy(h)
+hmask[mask_rho .== 0] .= NaN;
+pcolormesh(lon,lat,hmask);
+colorbar()
+## or colorbar(orientation="horizontal")
+gca().set_aspect(1/cosd(mean(lat)))
+
+title("smoothed bathymetry [m]");
+savefig("smoothed_bathymetry.png");
+
+# ![](smoothed_bathymetry.png)
+
+# ## Surface temperature
+
+# The surface surface temperature (or salinity) of the model output or
+# climatology file can be visualized as follows.
+# The parameter `n` is the time instance to plot.
+# Make sure that your current working directory
+# contains the file to plot (use e.g. `;cd ~/ROMS-implementation-test/` to plot `roms_his.nc`)
+
+
+## instance to plot
+n = 1
+
+ds = NCDataset("roms_his.nc")
+temp = nomissing(ds["temp"][:,:,end,n],NaN);
+temp[mask_rho .== 0] .= NaN;
+
+if haskey(ds,"time")
+    ## for the climatology file
+    time = ds["time"][:]
+else
+    ## for ROMS output
+    time = ds["ocean_time"][:]
+end
+
+clf();
+pcolormesh(lon,lat,temp)
+gca().set_aspect(1/cosd(mean(lat)))
+colorbar();
+title("sea surface temperature [°C]")
+savefig("SST.png");
+
+# ![](SST.png)
+
+# Exercise:
+# * Plot salinity
+# * Plot different time instance (`n`)
+# * Where do we specify that the surface values are to be plotted? Plot different layers.
+
+
+# ## Surface velocity and elevation
+
+zeta = nomissing(ds["zeta"][:,:,n],NaN)
+u = nomissing(ds["u"][:,:,end,n],NaN);
+v = nomissing(ds["v"][:,:,end,n],NaN);
+
+mask_u = ds_grid["mask_u"][:,:];
+mask_v = ds_grid["mask_v"][:,:];
+
+u[mask_u .== 0] .= NaN;
+v[mask_v .== 0] .= NaN;
+zeta[mask_rho .== 0] .= NaN;
+
+## ROMS uses an Arakawa C grid
+u_r = cat(u[1:1,:], (u[2:end,:] .+ u[1:end-1,:])/2, u[end:end,:], dims=1);
+v_r = cat(v[:,1:1], (v[:,2:end] .+ v[:,1:end-1])/2, v[:,end:end], dims=2);
+
+## all sizes should be the same
+size(u_r), size(v_r), size(mask_rho)
+
+clf();
+pcolormesh(lon,lat,zeta)
+colorbar();
+## plot only a single arrow for r x r grid cells
+r = 3;
+i = 1:r:size(lon,1);
+j = 1:r:size(lon,2);
+q = quiver(lon[i,j],lat[i,j],u_r[i,j],v_r[i,j])
+quiverkey(q,0.9,0.85,1,"1 m/s",coordinates="axes")
+title("surface currents [m/s] and elevation [m]");
+gca().set_aspect(1/cosd(mean(lat)))
+savefig("surface_zeta_uv.png");
+
+# ![](surface_zeta_uv.png)
+
+# Exercise:
+# * The surface currents seems to follow lines of constant surface elevation. Explain why this is to be expected.
+
+# ## Vertical section
+
+# In this example we will plot a vertical section by slicing the
+# model output at a given index.
+
+# It is very important that the parameters (`opt`) defining the vertical layer match the parameters values choosen when ROMS was setup.
+
+opt = (
+    Tcline = 50,   # m
+    theta_s = 5,   # surface refinement
+    theta_b = 0.4, # bottom refinement
+    nlevels = 32,  # number of vertical levels
+    Vtransform  = 2,
+    Vstretching = 4,
+)
+
+hmin = minimum(h)
+hc = min(hmin,opt.Tcline)
+z_r = ROMS.set_depth(opt.Vtransform, opt.Vstretching,
+                   opt.theta_s, opt.theta_b, hc, opt.nlevels,
+                   1, h);
+
+temp = nomissing(ds["temp"][:,:,:,n],NaN);
+
+mask3 = repeat(mask_rho,inner=(1,1,opt.nlevels))
+lon3 = repeat(lon,inner=(1,1,opt.nlevels))
+lat3 = repeat(lat,inner=(1,1,opt.nlevels))
+
+temp[mask3 .== 0] .= NaN;
+
+i = 20;
+
+clf()
+contourf(lat3[i,:,:],z_r[i,:,:],temp[i,:,:],40)
+ylim(-300,0);
+xlabel("latitude")
+ylabel("depth [m]")
+title("temperature at $(round(lon[i,1],sigdigits=4)) °E")
+colorbar();
+
+ax2 = gcf().add_axes([0.1,0.18,0.4,0.3]) # inset plot
+ax2.pcolormesh(lon,lat,temp[:,:,end])
+ax2.set_aspect(1/cosd(mean(lat)))
+ax2.plot(lon[i,[1,end]],lat[i,[1,end]],"m")
+
+savefig("temp_section1.png");
+
+# ![temp_section1](temp_section1.png)
+
+# Exercise:
+# * Plot a section at different longitude and latitude
+
+# ## Horizontal section
+
+# A horizontal at the fixed depth of 200 m is extracted and plotted.
+
+tempi = ROMS.model_interp3(lon,lat,z_r,temp,lon,lat,[-200])
+mlon,mlat,mdata = GeoDatasets.landseamask(resolution='f', grid=1.25)
+
+clf();
+pcolormesh(lon,lat,tempi[:,:,1])
+colorbar();
+ax = axis()
+contourf(mlon,mlat,mdata',[0.5, 3],colors=["gray"])
+axis(ax)
+gca().set_aspect(1/cosd(mean(lat)))
+title("temperature at 200 m [°C]")
+savefig("temp_hsection_200.png");
+
+# ![](temp_hsection_200.png)
+
+# ## Arbitrary vertical section
+
+# The vectors `section_lon` and `section_lat` define the coordinates where we want to extract
+# the surface temperature.
+
+
+section_lon = LinRange(8.18, 8.7,100);
+section_lat = LinRange(43.95, 43.53,100);
+
+using Interpolations
+
+function section_interp(v)
+    itp = interpolate((lon[:,1],lat[1,:]),v,Gridded(Linear()))
+    return itp.(section_lon,section_lat)
+end
+
+section_temp = mapslices(section_interp,temp,dims=(1,2))
+section_z = mapslices(section_interp,z_r,dims=(1,2))
+
+section_x = repeat(section_lon,inner=(1,size(temp,3)))
+
+clf()
+contourf(section_x,section_z[:,1,:],section_temp[:,1,:],50)
+ylim(-500,0)
+colorbar()
+xlabel("longitude")
+ylabel("depth")
+title("temperature section [°C]");
+
+ax2 = gcf().add_axes([0.4,0.2,0.4,0.3])
+ax2.pcolormesh(lon,lat,temp[:,:,end])
+axis("on")
+ax2.set_aspect(1/cosd(mean(lat)))
+ax2.plot(section_lon,section_lat,"m")
+
+savefig("temp_vsection.png");
+
+# ![](temp_vsection.png)