Merge pull request #72 from gaelforget/v0.2.28a

V0.2.28a
JuliaClimate · Nov 26, 2021 · 42fde35 · 42fde35
2 parents 0b134bf + 0f4ea41
commit 42fde35
Show file tree

Hide file tree

Showing 6 changed files with 217 additions and 89 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "MeshArrays"
 uuid = "cb8c808f-1acf-59a3-9d2b-6e38d009f683"
 authors = ["gaelforget <[email protected]>"]
-version = "0.2.27"
+version = "0.2.28"
 
 [deps]
 CatViews = "81a5f4ea-a946-549a-aa7e-2a7f63a27d31"

diff --git a/docs/src/API.md b/docs/src/API.md
@@ -15,9 +15,11 @@ varmeta
 
 ```@docs
 GridSpec
+GridLoad
+GridLoadVar
+UnitGrid
 GridOfOnes
 simple_periodic_domain
-GridLoad
 MeshArrays.read
 MeshArrays.write
 exchange

diff --git a/src/GridPaths.jl b/src/GridPaths.jl
@@ -0,0 +1,89 @@
+
+#bring end points to the equator -> define 3D rotation matrix
+function rotate_points(lons,lats)
+	#... and of end points
+	x0=cosd.(lats).*cosd.(lons)
+	y0=cosd.(lats).*sind.(lons)
+	z0=sind.(lats)
+
+	#get the rotation matrix:
+	#1) rotate around x axis to put first point at z=0
+	theta=atan(-z0[1],y0[1])
+	R1=[[1;0;0] [0;cos(theta);sin(theta)] [0;-sin(theta);cos(theta)]]
+	tmp0=[x0;y0;z0]; tmp1=R1*tmp0; x1=tmp1[1,:]; y1=tmp1[2,:]; z1=tmp1[3,:]
+	x0=x1; y0=y1; z0=z1
+	#2) rotate around z axis to put first point at y=0
+	theta=atan(x0[1],y0[1])
+	R2=[[cos(theta);sin(theta);0] [-sin(theta);cos(theta);0] [0;0;1]]
+	tmp0=transpose([x0 y0 z0]); tmp1=R2*tmp0; x1=tmp1[1,:]; y1=tmp1[2,:]; z1=tmp1[3,:]
+	x0=x1; y0=y1; z0=z1
+	#3) rotate around y axis to put second point at z=0
+	theta=atan(-z0[2],-x0[2])
+	R3=[[cos(theta);0;-sin(theta)] [0;1;0] [sin(theta);0;cos(theta)]]
+	tmp0=transpose([x0 y0 z0]); tmp1=R3*tmp0; x1=tmp1[1,:]; y1=tmp1[2,:]; z1=tmp1[3,:]
+	x0=x1; y0=y1; z0=z1
+
+	x0,y0,z0,R3*R2*R1
+end
+
+function rotate_XCYC(Γ,R)
+	#3D carthesian coordinates:
+	lon=Γ.XC; lat=Γ.YC
+	x=cosd.(lat)*cosd.(lon)
+	y=cosd.(lat)*sind.(lon)
+	z=sind.(lat)
+
+	#rotation R:
+    γ=Γ.XC.grid
+	tmpx=γ.write(x); tmpy=γ.write(y); tmpz=γ.write(z)
+	tmp1=findall((!isnan).(tmpx))
+	tmpx2=tmpx[tmp1]; tmpy2=tmpy[tmp1]; tmpz2=tmpz[tmp1]
+	tmp2=[tmpx2';tmpy2';tmpz2']
+
+	tmp3=R*tmp2
+
+	tmpx2=tmp3[1,:]; tmpy2=tmp3[2,:]; tmpz2=tmp3[3,:]
+	tmpx[tmp1]=tmpx2; tmpy[tmp1]=tmpy2; tmpz[tmp1]=tmpz2
+	x=γ.read(tmpx,Γ.XC); y=γ.read(tmpy,Γ.XC); z=γ.read(tmpz,Γ.XC)
+
+	x,y,z
+end
+
+##
+
+function shorter_paths!(xyz,xyz0,msk_in)
+	(x0,y0,z0)=xyz0[:]
+	(x,y,z)=xyz[:]
+
+	#split in two segments:
+	theta=zeros(2)
+	theta[1]=atan(y0[1],x0[1])
+	theta[2]=atan(y0[2],x0[2])
+
+    γ=msk_in[1].grid
+	tmpx=γ.write(x); tmpy=γ.write(y); tmpz=γ.write(z);
+	tmptheta=atan.(tmpy,tmpx)
+	if theta[2]<0;
+		tmp00=findall(tmptheta.<=theta[2])
+		tmptheta[tmp00].=tmptheta[tmp00]+2*pi
+		theta[2]=theta[2]+2*pi
+	end
+
+	msk_out=[]
+	for kk in 1:3
+		#select field to treat:
+		mm=msk_in[kk]
+		#select the shorther segment:
+		tmpm=γ.write(mm)
+       if theta[2]-theta[1]<=pi
+            tmpm[findall( (tmptheta.>theta[2]).|(tmptheta.<theta[1]) )].=0.0
+        else
+            tmpm[findall( (tmptheta.<=theta[2]).&(tmptheta.>=theta[1]) )].=0.0
+        end
+        mm=γ.read(tmpm,mm);
+        #store result:
+		push!(msk_out,mm)
+    end
+	msk_out[:]
+end
+
diff --git a/src/Grids.jl b/src/Grids.jl
@@ -212,72 +212,82 @@ function GridLoad(γ::gcmgrid;option="minimal")
 
     Γ=Dict()
 
-    pc=fill(0.5,2); pg=fill(0.0,2); pu=[0.,0.5]; pv=[0.5,0.];
     if option=="full"
         list_n=("XC","XG","YC","YG","RAC","RAW","RAS","RAZ","DXC","DXG","DYC","DYG","Depth");
-        list_u=(u"°",u"°",u"°",u"°",u"m^2",u"m^2",u"m^2",u"m^2",u"m",u"m",u"m",u"m",u"m")
-        list_p=(pc,pg,pc,pg,pc,pu,pv,pg,pu,pv,pv,pu,pc)
+        if (!isempty(filter(x -> occursin("AngleCS",x), readdir(γ.path))))
+            list_n=(list_n...,"AngleCS","AngleSN");
+        end
+        list_n=(list_n...,"DRC","DRF","RC","RF")
+        list_n=(list_n...,"hFacC","hFacS","hFacW");
     else
         list_n=("XC","YC");
-        list_u=(u"°",u"°")
-        list_p=(pc,pc)
     end
 
-    if (option=="full")&&(!isempty(filter(x -> occursin("AngleCS",x), readdir(γ.path))))
-        list_n=(list_n...,"AngleCS","AngleSN");
-        list_u=(list_u...,1.0,1.0)
-        list_p=(list_p...,pc,pc)
-    end
+    [Γ[ii]=GridLoadVar(ii,γ) for ii in list_n]
+    option=="full" ? GridAddWS!(Γ) : nothing
+    return Dict_to_NamedTuple(Γ)
+end
 
-    for ii=1:length(list_n)
-        m=varmeta(list_u[ii],list_p[ii],missing,list_n[ii],list_n[ii])
-        tmp1=γ.read(joinpath(γ.path,list_n[ii]*".data"),MeshArray(γ,γ.ioPrec;meta=m))
-        tmp2=Symbol(list_n[ii])
-        @eval (($tmp2) = ($tmp1))
-        Γ[list_n[ii]]=tmp1
-    end
+"""
+    GridLoadVar(nam::String,γ::gcmgrid)
 
-    γ.ioPrec==Float64 ? reclen=8 : reclen=4
+Return a grid variable read from files located in `γ.path` (see `?GridSpec`, `?GridLoad`).
 
-    if option=="full"
-        list_n=("DRC","DRF","RC","RF")
-    else
-        list_n=()
-    end
-    for ii=1:length(list_n)
-        fil=joinpath(γ.path,list_n[ii]*".data")
-        tmp1=stat(fil)
-        n3=Int64(tmp1.size/reclen)
+Based on the MITgcm naming convention, grid variables are:
 
-        fid = open(fil)
-        tmp1 = Array{γ.ioPrec,1}(undef,n3)
-        read!(fid,tmp1)
-        tmp1 = hton.(tmp1)
+- XC, XG, YC, YG, AngleCS, AngleSN, hFacC, hFacS, hFacW, Depth.
+- RAC, RAW, RAS, RAZ, DXC, DXG, DYC, DYG.
+- DRC, DRF, RC, RF (one-dimensional)
 
-        tmp2=Symbol(list_n[ii])
-        @eval (($tmp2) = ($tmp1))
-        Γ[list_n[ii]]=tmp1
-    end
+```jldoctest
+using MeshArrays
 
-    f=readdir(γ.path)
-    if (option=="full")&&(sum(occursin.("hFacC",f))>0)
-        list_n=("hFacC","hFacS","hFacW");
-        list_u=(1.0,1.0,1.0)
-        list_p=(fill(0.5,3),[0.,0.5,0.5],[0.5,0.,0.5])
-        n3=length(Γ["RC"])
-        for ii=1:length(list_n)
-            m=varmeta(list_u[ii],list_p[ii],missing,list_n[ii],list_n[ii]);
-            tmp1=γ.read(joinpath(γ.path,list_n[ii]*".data"),MeshArray(γ,γ.ioPrec,n3;meta=m))
-            tmp2=Symbol(list_n[ii])
-            @eval (($tmp2) = ($tmp1))
-            Γ[list_n[ii]]=tmp1
-        end
-    end
+γ = GridSpec("CubeSphere",MeshArrays.GRID_CS32)
+XC = GridLoadVar("XC",γ)
 
-    option=="full" ? GridAddWS!(Γ) : nothing
+isa(XC,MeshArray)
 
-    return Dict_to_NamedTuple(Γ)
+# output
 
+true
+```
+"""
+function GridLoadVar(nam::String,γ::gcmgrid)
+    pc=fill(0.5,2); pg=fill(0.0,2); pu=[0.,0.5]; pv=[0.5,0.];
+    list_n=("XC","XG","YC","YG","RAC","RAW","RAS","RAZ","DXC","DXG","DYC","DYG","Depth","AngleCS","AngleSN")
+    list_u=(u"°",u"°",u"°",u"°",u"m^2",u"m^2",u"m^2",u"m^2",u"m",u"m",u"m",u"m",u"m",1.0,1.0)
+    list_p=(pc,pg,pc,pg,pc,pu,pv,pg,pu,pv,pv,pu,pc,pc,pc)
+    #
+    list3d_n=("hFacC","hFacS","hFacW");
+    list3d_u=(1.0,1.0,1.0)
+    list3d_p=(fill(0.5,3),[0.,0.5,0.5],[0.5,0.,0.5])
+    #
+    list1d_n=("DRC","DRF","RC","RF")
+
+    if sum(nam.==list_n)==1
+        ii=findall(nam.==list_n)[1]
+        m=varmeta(list_u[ii],list_p[ii],missing,list_n[ii],list_n[ii])
+        tmp1=γ.read(joinpath(γ.path,list_n[ii]*".data"),MeshArray(γ,γ.ioPrec;meta=m))
+    elseif sum(nam.==list1d_n)==1
+        fil=joinpath(γ.path,nam*".data")
+        γ.ioPrec==Float64 ? reclen=8 : reclen=4
+        n3=Int64(stat(fil).size/reclen)
+
+        fid = open(fil)
+        tmp1 = Array{γ.ioPrec,1}(undef,n3)
+        read!(fid,tmp1)
+        tmp1 = hton.(tmp1)
+    elseif sum(nam.==list3d_n)==1
+        fil=joinpath(γ.path,"RC.data")
+        γ.ioPrec==Float64 ? reclen=8 : reclen=4
+        n3=Int64(stat(fil).size/reclen)
+
+        ii=findall(nam.==list3d_n)[1]
+        m=varmeta(list3d_u[ii],list3d_p[ii],missing,list3d_n[ii],list3d_n[ii]);
+        tmp1=γ.read(joinpath(γ.path,list3d_n[ii]*".data"),MeshArray(γ,γ.ioPrec,n3;meta=m))
+    else
+        tmp1=missing
+    end
 end
 
 """

diff --git a/src/MeshArrays.jl b/src/MeshArrays.jl
@@ -22,6 +22,7 @@ GRID_CS32_download() = artifact"GRID_CS32"
 
 include("Types.jl")
 include("Grids.jl")
+include("GridPaths.jl")
 include("Operations.jl")
 include("Exchanges.jl")
 include("ReadWrite.jl")
@@ -32,11 +33,12 @@ include("VerticalDimension.jl")
 
 export AbstractMeshArray, MeshArray, InnerArray, OuterArray, varmeta
 export gcmgrid, exchange, convergence, smooth, mask, isosurface
-export UnitGrid, simple_periodic_domain, GridSpec, GridLoad, GridOfOnes, GridAddWS!
+export UnitGrid, simple_periodic_domain, GridOfOnes 
+export GridSpec, GridLoad, GridLoadVar, GridAddWS!
 export Tiles, Interpolate, InterpolationFactors, knn
 export ScalarPotential, VectorPotential, ThroughFlow
 export UVtoUEVN, UVtoTransport, gradient, curl
-export StereographicProjection, LatitudeCircles
+export StereographicProjection, LatitudeCircles, Transect
 
 export location_is_out, NeighborTileIndices_dpdo, NeighborTileIndices_cs, RelocationFunctions_cs
 export update_location_cs!, update_location_llc!, update_location_dpdo!

diff --git a/src/Operations.jl b/src/Operations.jl
@@ -313,48 +313,73 @@ end
 Compute integration paths that follow latitude circles
 """
 function LatitudeCircles(LatValues,Γ::NamedTuple)
-
     LatitudeCircles=Array{NamedTuple}(undef,length(LatValues))
-
     for j=1:length(LatValues)
         mskCint=1*(Γ.YC .>= LatValues[j])
-        mskC=similar(mskCint)
-        mskW=similar(mskCint)
-        mskS=similar(mskCint)
-
-        mskCint=exchange(mskCint,1)
-
-        for i=1:mskCint.grid.nFaces
-            tmp1=mskCint.f[i]
-            # tracer mask:
-            tmp2=tmp1[2:end-1,1:end-2]+tmp1[2:end-1,3:end]+
-            tmp1[1:end-2,2:end-1]+tmp1[3:end,2:end-1]
-            mskC.f[i]=1((tmp2.>0).&(tmp1[2:end-1,2:end-1].==0))
-            # velocity masks:
-            mskW.f[i]=tmp1[2:end-1,2:end-1] - tmp1[1:end-2,2:end-1]
-            mskS.f[i]=tmp1[2:end-1,2:end-1] - tmp1[2:end-1,1:end-2]
-        end
-
-        function MskToTab(msk::MeshArray)
-          n=Int(sum(msk .!= 0)); k=0
-          tab=Array{Int,2}(undef,n,4)
-          for i=1:msk.grid.nFaces
-            a=msk.f[i]
-            b=findall( a .!= 0)
-            for ii in eachindex(b)
-              k += 1
-              tab[k,:]=[i,b[ii][1],b[ii][2],a[b[ii]]]
-            end
-          end
-          return tab
-        end
-
+        mskC,mskW,mskS=edge_mask(mskCint)
         LatitudeCircles[j]=(lat=LatValues[j],
         tabC=MskToTab(mskC),tabW=MskToTab(mskW),tabS=MskToTab(mskS))
     end
-
     return LatitudeCircles
+end
+
+function MskToTab(msk::MeshArray)
+  n=Int(sum(msk .!= 0)); k=0
+  tab=Array{Int,2}(undef,n,4)
+  for i in eachindex(msk)
+    a=msk[i]
+    b=findall( a .!= 0)
+    for ii in eachindex(b)
+      k += 1
+      tab[k,:]=[i,b[ii][1],b[ii][2],a[b[ii]]]
+    end
+  end
+  return tab
+end
+
+function edge_mask(mskCint::MeshArray)
+  mskCint=1.0*mskCint
+
+  #treat the case of blank tiles:
+  #mskCint[findall(RAC.==0)].=NaN
+
+  mskC=similar(mskCint)
+  mskW=similar(mskCint)
+  mskS=similar(mskCint)
+
+  mskCint=exchange(mskCint,1)
+
+  for i in eachindex(mskCint)
+      tmp1=mskCint[i]
+      # tracer mask:
+      tmp2=tmp1[2:end-1,1:end-2]+tmp1[2:end-1,3:end]+
+      tmp1[1:end-2,2:end-1]+tmp1[3:end,2:end-1]
+      mskC[i]=1((tmp2.>0).&(tmp1[2:end-1,2:end-1].==0))
+      # velocity masks:
+      mskW[i]=tmp1[2:end-1,2:end-1] - tmp1[1:end-2,2:end-1]
+      mskS[i]=tmp1[2:end-1,2:end-1] - tmp1[2:end-1,1:end-2]
+  end
+
+  #treat the case of blank tiles:
+  #mskC[findall(isnan.(mskC))].=0.0
+  #mskW[findall(isnan.(mskW))].=0.0
+  #mskS[findall(isnan.(mskS))].=0.0
+
+  return mskC,mskW,mskS
+end
+
+##
 
+function Transect(name,lons,lats,Γ)
+  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))
+  tabC=MskToTab(mskCedge)
+  tabW=MskToTab(mskWedge)
+  tabS=MskToTab(mskSedge)
+  return (name=name,tabC=tabC,tabW=tabW,tabS=tabS)
 end
 
 ##