Merge pull request #143 from gaelforget/v0p3p12_e

methods improvements

Author: gaelforget

ext/MeshArraysMakieExt.jl

@@ -739,6 +739,14 @@ function scatter!(pr_ax::PrAxis,lon,lat,kargs...; kwargs...)
 	scatter!(pr_ax.ax,x,y,kargs...; kwargs...)
 end
 
+plot(x::MeshArray; kwargs...) = begin
+	if ndims(x) == 1
+		heatmap(x; kwargs...)
+	else
+		heatmap(x[:,1]; kwargs...)
+	end
+end
+
 ##
 
 function plot(x::Union{gridpath,Vector{gridpath}})

src/Grids.jl

@@ -130,27 +130,39 @@ function simple_periodic_domain(np::Integer,nq=missing)
     return UnitGrid(γ)
 end
 
-## GridSpec function with default GridName argument:
+## GridSpec function with category argument:
 
-GridSpec() = GridSpec("PeriodicDomain","./")
+"""
+    GridSpec(category="PeriodicDomain",path=tempname(); ID=:unknown)
 
-## GridSpec function with GridName argument:
+- Select one of the pre-defined grids either by ID (keyword) or by category.
+- Return the corresponding `gcmgrid` specification, including the path where grid files can be accessed (`path`).
 
-"""
-    GridSpec(GridName,GridParentDir="./")
+1. selection by `ID`
 
-Select one of the pre-defined grids (by `GridName`) and return 
-the corresponding `gcmgrid` -- a global grid specification 
-which contains the grid files location (`GridParentDir`).
-    
+- `:LLC90`
+- `:CS32`
+- `:onedegree`
+- `:default`
+
+Example:
+
+```
+using MeshArrays
+g = GridSpec(ID=:LLC90)
+```
+
+note : the path to these fully supported grids are handled internally in `MeshArrays.jl`.
 
-Possible choices for `GridName`:
+2. by `category` and `path`
 
 - `"PeriodicDomain"`
 - `"PeriodicChannel"`
 - `"CubeSphere"`
 - `"LatLonCap"``
 
+Examples:
+
 ```jldoctest; output = false
 using MeshArrays
 g = GridSpec()
@@ -164,51 +176,62 @@ isa(g,gcmgrid)
 true
 ```
 """
-function GridSpec(GridName,GridParentDir="./")
+function GridSpec(category="PeriodicDomain",path=tempname(); ID=:unknown)
 
-grDir=GridParentDir
-if GridName=="LatLonCap"
+if category=="LatLonCap"
     nFaces=5
     grTopo="LatLonCap"
     ioSize=[90 1170]
     facesSize=[(90, 270), (90, 270), (90, 90), (270, 90), (270, 90)]
     ioPrec=Float64
-elseif GridName=="CubeSphere"
+elseif category=="CubeSphere"
     nFaces=6
     grTopo="CubeSphere"
     ioSize=[32 192]
     facesSize=[(32, 32), (32, 32), (32, 32), (32, 32), (32, 32), (32, 32)]
     ioPrec=Float32
-elseif GridName=="PeriodicChannel"
+elseif category=="PeriodicChannel"
     nFaces=1
     grTopo="PeriodicChannel"
     ioSize=[360 160]
     facesSize=[(360, 160)]
     ioPrec=Float32
-elseif GridName=="PeriodicDomain"
+elseif category=="PeriodicDomain"
     nFaces=4
     grTopo="PeriodicDomain"
     ioSize=[80 42]
     facesSize=[(40, 21), (40, 21), (40, 21), (40, 21)]
     ioPrec=Float32
 else
-    error("unknown GridName case")
+    error("unknown category case")
 end
 
-return gcmgrid(grDir,grTopo,nFaces,facesSize, ioSize, ioPrec, read, write)
+if ID==:unknown
+    gcmgrid(path,grTopo,nFaces,facesSize, ioSize, ioPrec, read, write)
+elseif ID==:LLC90
+    GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
+elseif ID==:CS32
+    GridSpec("CubeSphere",MeshArrays.GRID_CS32)
+elseif ID==:onedegree
+    GridSpec("PeriodicChannel",MeshArrays.GRID_LL360)
+elseif ID==:default
+    GridSpec()
+else
+    error("unknwown grid")
+end
 
 end
 
 ## GridLoad function
 
 """
-    GridLoad(γ::gcmgrid;option="minimal")
+    GridLoad(γ=GridSpec(); ID=:default, option=:minimal)
 
-Return a `NamedTuple` of grid variables read from files located in `γ.path` (see `?GridSpec`).
-
-By default, option="minimal" means that only grid cell center positions (XC, YC) are loaded. 
-
-The "full" option provides a complete set of grid variables. 
+- Return a `NamedTuple` of grid variables read from files located in `γ.path` (see `?GridSpec`).
+- option : 
+  - (default) option=:minimal means that only grid cell center positions (XC, YC) are loaded. 
+  - option=:full provides a complete set of 2D grid variables. 
+  - option=:full provides a complete set of 2D & 3d grid variables. 
 
 Based on the MITgcm naming convention, grid variables are:
 
@@ -222,11 +245,7 @@ https://mitgcm.readthedocs.io/en/latest/algorithm/algorithm.html#spatial-discret
 
 ```jldoctest; output = false
 using MeshArrays
-
-γ = GridSpec("CubeSphere",MeshArrays.GRID_CS32)
-#γ = GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
-#γ = GridSpec("PeriodicChannel",MeshArrays.GRID_LL360)
-
+γ = GridSpec(ID=:LLC90)
 Γ = GridLoad(γ;option="full")
 
 isa(Γ.XC,MeshArray)
@@ -236,33 +255,38 @@ isa(Γ.XC,MeshArray)
 true
 ```
 """
-function GridLoad(γ::gcmgrid;option="minimal")
+function GridLoad(γ=GridSpec(); ID=:default, option=:minimal)
+
+    gr = (ID!==:default ? GridSpec(ID=ID) : γ)
 
-    γ.path==GRID_CS32 ? GRID_CS32_download() : nothing
-    γ.path==GRID_LL360 ? GRID_LL360_download() : nothing
-    γ.path==GRID_LLC90 ? GRID_LLC90_download() : nothing
+    gr.path==GRID_CS32 ? GRID_CS32_download() : nothing
+    gr.path==GRID_LL360 ? GRID_LL360_download() : nothing
+    gr.path==GRID_LLC90 ? GRID_LLC90_download() : nothing
 
     Γ=Dict()
 
-    if option=="full"
-        list_n=("XC","XG","YC","YG","RAC","RAW","RAS","RAZ","DXC","DXG","DYC","DYG","Depth");
-        if (!isempty(filter(x -> occursin("AngleCS",x), readdir(γ.path))))
+    op=string(option)
+    if op=="full"
+        list_n=("XC","XG","YC","YG","RAC","RAW","RAS","RAZ","DXC","DXG","DYC","DYG","Depth")
+        if (!isempty(filter(x -> occursin("AngleCS",x), readdir(gr.path))))
             list_n=(list_n...,"AngleCS","AngleSN");
         end
         list_n=(list_n...,"DRC","DRF","RC","RF")
-        list_n=(list_n...,"hFacC","hFacS","hFacW");
-    elseif option=="light"
-        list_n=("XC","XG","YC","YG","RAC","DXC","DXG","DYC","DYG","Depth");
-        if (!isempty(filter(x -> occursin("AngleCS",x), readdir(γ.path))))
-            list_n=(list_n...,"AngleCS","AngleSN");
+        list_n=(list_n...,"hFacC","hFacS","hFacW")
+    elseif op=="light"
+        list_n=("XC","XG","YC","YG","RAC","DXC","DXG","DYC","DYG","Depth")
+        if (!isempty(filter(x -> occursin("AngleCS",x), readdir(gr.path))))
+            list_n=(list_n...,"AngleCS","AngleSN")
         end
         list_n=(list_n...,"DRC","DRF","RC","RF")
+    elseif op=="minimal"||op=="minimum"
+        list_n=("XC","YC")
     else
-        list_n=("XC","YC");
+        error("unknown option")
     end
 
-    [Γ[ii]=GridLoadVar(ii,γ) for ii in list_n]
-    option=="full" ? GridAddWS!(Γ) : nothing
+    [Γ[ii]=GridLoadVar(ii,gr) for ii in list_n]
+    option=="full"||option=="light" ? GridAddWS!(Γ) : nothing
     return Dict_to_NamedTuple(Γ)
 end
 

src/Operations.jl

@@ -331,26 +331,45 @@ end
 ## LatitudeCircles function
 
 """
-    LatitudeCircles(LatValues,Γ::NamedTuple; format=:gridpath)
+    LatitudeCircles(LatValues,Γ::NamedTuple; format=:gridpath, range=(0.0,360.0))
 
-Compute integration paths that follow latitude circles
+Compute integration paths that follow latitude circles, within the specified longitude `range`.
 """
-function LatitudeCircles(LatValues,Γ::NamedTuple; format=:gridpath)
+function LatitudeCircles(LatValues,Γ::NamedTuple; 
+  format=:gridpath, range=(0.0,360.0))
   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]=
-          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
+        LatitudeCircles[j]=LatitudeCircle(LatValues[j],Γ; 
+        format=format,range=range)
     end
-    return LatitudeCircles
+    (length(LatValues)==1 ? LatitudeCircles[1] : LatitudeCircles)
+end
+
+function LatitudeCircle(lat,Γ::NamedTuple; 
+  format=:gridpath, range=(0.0,360.0))
+      mskCint=1*(Γ.YC .>= lat)
+      mskC,mskW,mskS=edge_mask(mskCint)
+      restrict_longitudes!(mskC,Γ.XC,range=range)
+      restrict_longitudes!(mskS,Γ.XS,range=range)
+      restrict_longitudes!(mskW,Γ.XW,range=range)
+      LC=if format==:NamedTuple
+        (lat=LatValues[j],tabC=MskToTab(mskC),
+        tabW=MskToTab(mskW),tabS=MskToTab(mskS))
+      else
+        gridpath(name="Parallel $lat", grid=Γ,
+        C=MskToTab(mskC),W=MskToTab(mskW),S=MskToTab(mskS))
+      end
+
+end
+
+is_in_lon_range(x,range)=(range[2].-range[1]>=360)||
+  (mod(x-range[1],360).<mod(range[2].-range[1],360))
+
+function restrict_longitudes!(x::MeshArray,lon::MeshArray;range=(0.0,360.0))
+  for f in 1:x.grid.nFaces
+    x[f].=x[f].*is_in_lon_range.(lon[f],Ref(range))
+  end
 end
 
 function MskToTab(msk::MeshArray)