Merge pull request #142 from gaelforget/v0p3p13b

V0p3p13b

Project.toml

@@ -1,7 +1,7 @@
 name = "MeshArrays"
 uuid = "cb8c808f-1acf-59a3-9d2b-6e38d009f683"
 authors = ["gaelforget <gforget@mit.edu>"]
-version = "0.3.11"
+version = "0.3.12"
 
 [deps]
 CatViews = "81a5f4ea-a946-549a-aa7e-2a7f63a27d31"

ext/MeshArraysMakieExt.jl

@@ -9,7 +9,11 @@ import MeshArrays: plot_examples
 import MeshArrays: ProjAxis
 import MeshArrays: grid_lines!
 
-import Makie: heatmap, scatter, scatter!, surface!, lines!, heatmap!, contour!, contourf!
+import MeshArrays: gridpath
+
+import Makie: plot, plot!, heatmap, scatter, scatter!, surface!
+import Makie: lines!, heatmap!, contour!, contourf!
+
 LineString=Makie.LineString
 Observable=Makie.Observable
 GeometryBasics=Makie.GeometryBasics
@@ -735,5 +739,31 @@ function scatter!(pr_ax::PrAxis,lon,lat,kargs...; kwargs...)
 	scatter!(pr_ax.ax,x,y,kargs...; kwargs...)
 end
 
+##
+
+function plot(x::Union{gridpath,Vector{gridpath}})
+	fig=Figure(); ax=Axis(fig[1,1],limits=(-180.0,180.0,-90.0,90.0))
+	if isa(x,gridpath)
+		plot!(x)
+	else
+		[plot!(y) for y in x]
+	end
+	fig
+end
+
+function plot!(x::gridpath)
+	np=size(x.C,1)
+	lon=zeros(np)
+	lat=zeros(np)
+	for p in 1:np
+		f=x.C[p,1]; i=x.C[p,2]; j=x.C[p,3]; q=x.C[p,4]
+		lon[p]=x.grid.XC[f][i,j]
+		lat[p]=x.grid.YC[f][i,j]
+	end
+	scatter!(lon,lat)
+end
+
 end # module
 
+
+

src/MeshArrays.jl

@@ -26,7 +26,7 @@ include("Interpolation.jl")
 include("VerticalDimension.jl")
 include("demo.jl")
 
-export AbstractMeshArray, MeshArray, gcmgrid, varmeta
+export AbstractMeshArray, MeshArray, gcmgrid, varmeta, gridpath
 export GridSpec, GridLoad, GridLoadVar, UnitGrid, simple_periodic_domain
 export exchange, Tiles, Tiles!, Interpolate, InterpolationFactors, knn, interpolation_setup
 #The following exch_UV differs from normal exchange; incl. exch_UV_N

src/Operations.jl

@@ -296,8 +296,8 @@ function ThroughFlow(VectorField,IntegralPath,Γ::NamedTuple)
         #do_dz==1 ? mskS=Γ.DRF[i3]*mskS : nothing
         #
         #method 2: less slow
-        tabW=IntegralPath.tabW
-        tabS=IntegralPath.tabS
+        tabW=(isa(IntegralPath,NamedTuple) ? IntegralPath.tabW : IntegralPath.W)
+        tabS=(isa(IntegralPath,NamedTuple) ? IntegralPath.tabS : IntegralPath.S)
         for i4=1:n[4]
             #method 1: quite slow
             #trsp[1,i3,i4]=sum(mskW*U[:,:,i3,i4])+sum(mskS*V[:,:,i3,i4])
@@ -331,17 +331,24 @@ end
 ## LatitudeCircles function
 
 """
-    LatitudeCircles(LatValues,Γ::NamedTuple)
+    LatitudeCircles(LatValues,Γ::NamedTuple; format=:gridpath)
 
 Compute integration paths that follow latitude circles
 """
-function LatitudeCircles(LatValues,Γ::NamedTuple)
-    LatitudeCircles=Array{NamedTuple}(undef,length(LatValues))
+function LatitudeCircles(LatValues,Γ::NamedTuple; format=:gridpath)
+  T=(format==:NamedTuple ? NamedTuple : gridpath)
+  LatitudeCircles=Array{T}(undef,length(LatValues))
     for j=1:length(LatValues)
         mskCint=1*(Γ.YC .>= LatValues[j])
         mskC,mskW,mskS=edge_mask(mskCint)
-        LatitudeCircles[j]=(lat=LatValues[j],
-        tabC=MskToTab(mskC),tabW=MskToTab(mskW),tabS=MskToTab(mskS))
+        LatitudeCircles[j]=
+          if format==:NamedTuple
+            (lat=LatValues[j],tabC=MskToTab(mskC),
+            tabW=MskToTab(mskW),tabS=MskToTab(mskS))
+          else
+            gridpath(name="Parallel $(LatValues[j])", grid=Γ,
+            C=MskToTab(mskC),W=MskToTab(mskW),S=MskToTab(mskS))
+          end
     end
     return LatitudeCircles
 end
@@ -400,20 +407,36 @@ end
 ##
 
 """
-    Transect(name,lons,lats,Γ)
+    Transect(name,lons,lats,Γ; segment=:short, format=:gridpath)
 
 Compute integration paths that follow a great circle between two geolocations give by `lons`, `lats`.
+
 """
-function Transect(name,lons,lats,Γ)
+function Transect(name,lons,lats,Γ; segment=:short, format=:gridpath)
   x0,y0,z0,R=rotate_points(lons,lats)
   x,y,z=rotate_XCYC(Γ,R)
   mskCint=1.0*(z.>0)
   mskCedge,mskWedge,mskSedge=edge_mask(mskCint)
-  mskCedge,mskWedge,mskSedge=shorter_paths!((x,y,z),(x0,y0,z0),(mskCedge,mskWedge,mskSedge))
+
+  mskCedge,mskWedge,mskSedge=
+  if segment==:short
+    shorter_paths!((x,y,z),(x0,y0,z0),(mskCedge,mskWedge,mskSedge))
+  elseif segment==:long
+    C,W,S=shorter_paths!((x,y,z),(x0,y0,z0),(mskCedge,mskWedge,mskSedge))
+    C-mskCedge,W-mskWedge,S-mskSedge
+  elseif segment==:circle
+    mskCedge,mskWedge,mskSedge
+  end
+
   tabC=MskToTab(mskCedge)
   tabW=MskToTab(mskWedge)
   tabS=MskToTab(mskSedge)
-  return (name=name,tabC=tabC,tabW=tabW,tabS=tabS)
+
+  if format==:NamedTuple
+    (name=name,tabC=tabC,tabW=tabW,tabS=tabS)
+  else
+    gridpath(name=name,grid=Γ,C=tabC,W=tabW,S=tabS)
+  end
 end
 
 ##

src/Types.jl

@@ -209,3 +209,28 @@ function ones(a::AbstractMeshArray)
   [b.f[c].=1.0 for c in eachindex(b.f)]
   b
 end
+
+## derivative types
+
+
+"""
+    gridpath
+
+gridpath data structure.
+
+```
+using MeshArrays
+γ=GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
+Γ=GridLoad(γ;option="light")
+lons=[-68 -63]; lats=[-54 -66]; name="Drake Passage"
+Trsct=Transect(name,lons,lats,Γ)
+```
+"""
+Base.@kwdef struct gridpath
+#  grid::gcmgrid
+  grid::NamedTuple
+  name::String
+  C::Matrix
+  W::Matrix
+  S::Matrix
+end

test/runtests.jl

@@ -74,7 +74,7 @@ end
 
     #Meridional transport integral
     uv=Dict("U"=>Tx,"V"=>Ty,"dimensions"=>["x","y"])
-    L=-85.0:5.0:85.0; LC=LatitudeCircles(L,Γ)
+    L=-85.0:5.0:85.0; LC=LatitudeCircles(L,Γ,format=:gridpath)
     T=Array{Float64,1}(undef,length(LC))
     [T[i]=1e-6*ThroughFlow(uv,LC[i],Γ) for i=1:length(LC)]