Merge pull request #152 from gaelforget/v0p3p15a

V0p3p15a

Project.toml

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

docs/src/API.md

@@ -9,6 +9,8 @@ AbstractMeshArray
 MeshArrays.gcmarray
 gcmgrid
 varmeta
+gridpath
+gridmask
 ```
 
 More : 
@@ -21,7 +23,6 @@ MeshArrays.gcmfaces
 ## 2. Grids And I/O
 
 ```@docs
-UnitGrid
 simple_periodic_domain
 GridSpec
 GridLoad
@@ -37,6 +38,7 @@ exchange
 MeshArrays.read
 MeshArrays.read!
 MeshArrays.write
+UnitGrid
 ```
 
 ## 3. Interpolation
@@ -62,11 +64,19 @@ UVtoTransport
 UVtoUEVN
 ```
 
-## 5. Other
+## 5. Integration
+
+```@docs
+Integration.loops
+```
+
+## 6. Other
 
 ```@docs
 LatitudeCircles
 Transect
+demo.ocean_basins
+demo.extended_basin
 isosurface
-MA_datadep
+MeshArrays.mydatadep
 ```

ext/MeshArraysDataDepsExt.jl

@@ -2,7 +2,7 @@
 module MeshArraysDataDepsExt
 
    using DataDeps, MeshArrays
-   import MeshArrays: MA_datadep
+   import MeshArrays: mydatadep
 
 
    __init__() = begin
@@ -23,7 +23,7 @@ module MeshArraysDataDepsExt
 
 
    """
-       MA_datadep(nam="countries_shp1")
+       mydatadep(nam="countries_shp1")
        
    Download data dependency with predefined name; currently :
 
@@ -34,7 +34,7 @@ module MeshArraysDataDepsExt
    "interp_halfdeg"
    ```    
    """
-   MA_datadep(nam="countries_shp1") = begin
+   mydatadep(nam="countries_shp1") = begin
       withenv("DATADEPS_ALWAYS_ACCEPT"=>true) do
          if nam=="countries_shp1"
             datadep"countries_shp1"

src/Grids.jl

@@ -231,7 +231,7 @@ end
 - 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. 
+  - option=:full provides a complete set of 2D & 3D grid variables. 
 
 Based on the MITgcm naming convention, grid variables are:
 

src/Integration.jl

@@ -8,14 +8,14 @@ import MeshArrays: demo, MeshArray, gridmask
 
 ##
 
-example(;option=:hv,regions=:dlat_10,depths=[(0,7000)]) = begin
+example(;option=:loops,regions=:dlat_10,depths=[(0,7000)]) = begin
   g=GridSpec(ID=:LLC90); Γ=GridLoad(g)
   G=(hFacC=GridLoadVar("hFacC",g),RF=GridLoadVar("RF",g),
      RC=GridLoadVar("RC",g),RAC=GridLoadVar("RAC",g),
      DRF=GridLoadVar("DRF",g))
   G=merge(Γ,G)
 
-  M=define_boxes(option=option, regions=regions, grid=G, depths=depths)
+  M=define_sums(option=option, regions=regions, grid=G, depths=depths)
 
   diags=joinpath(pwd(),"diags")
   files=glob("state_3d_set1*.data",diags)
@@ -89,11 +89,11 @@ layer_mask(dF,d0,d1)=begin
 end
 
 """
-    define_boxes(;option=:hv,grid::NamedTuple, regions=:basins, depths=[(0,7000)])
+    define_sums(;option=:loops, grid::NamedTuple, regions=:basins, depths=[(0,7000)])
 
 Define regional integration function for each basin and depth range.
 """
-function define_boxes(;option=:hv, grid::NamedTuple, regions=:basins, depths=[(0,7000)])
+function define_sums(;option=:loops, grid::NamedTuple, regions=:basins, depths=[(0,7000)])
   dep=(isa(depths,Tuple) ? [depths] : depths)
   nd=length(dep)
   rgns=define_regions(option=regions,grid=grid) 
@@ -118,7 +118,7 @@ function define_boxes(;option=:hv, grid::NamedTuple, regions=:basins, depths=[(0
     tmp2d
   end
 
-  if option==:full
+  if option==:streamlined_loop
   #ocn_surf=[sum(xymsk(b)*(grid.hFacC[:,1].>0)*grid.RAC) for b in 1:nb]
   BX=(name=String[],volsum=Function[],volume=Float64[],
        ocn_surf=Float64[],tmp2d=tmp2d,tmp3d=tmp3d)
@@ -153,9 +153,9 @@ function define_boxes(;option=:hv, grid::NamedTuple, regions=:basins, depths=[(0
     push!(BXv.vint,f)
   end
 
-  if option==:hv
+  if option==:loops
     gridmask(rgns.map,BXh.name,depths,BXh.hsum,BXv.vint,tmp2d,tmp3d)
-  elseif option==:full
+  elseif option==:streamlined_loop
     gridmask(rgns.map,BX.name,depths,BX.volsum,[],tmp2d,tmp3d)
   else
     error("unknown option")
@@ -165,17 +165,18 @@ end
 #nonan(x)=[(isnan(y) ? 0.0 : y) for y in x]
 
 """
-    loop(boxes::NamedTuple; files=String[], var=:THETA, rd=read)
+    loops(mask::gridmask; files=String[], var=:THETA, rd=read)
 
 ```
 begin
-@everywhere using MeshArrays, MITgcm
-@everywhere rd(F,var,tmp)=read(read_mdsio(F,var),tmp)
-@everywhere G,M,files=Integration.example(option=:hv)
-	#,regions=(30,10),depths=Integration.DEPTHS)
+  @everywhere using MeshArrays, MITgcm
+  @everywhere rd(F,var,tim,tmp)=read(read_mdsio(F,var),tmp)
+  @everywhere G,M,files=Integration.example()
+    #,regions=(30,10),depths=Integration.DEPTHS)
 end;
-#H_3d=Integration.loop_3d(M,files=files,rd=rd)
-H_hv=Integration.loop_hv(M,files=files,rd=rd)
+
+H=Integration.loops(M,files=files,rd=rd)
+# Hbis=Integration.streamlined_loop(M,files=files,rd=rd)
 ```
 
 and to save results:
@@ -188,55 +189,59 @@ jldsave(output_path; depths=M.depths, integral=H, volume=vol, name=M.names)
 where vol is calculated as follows:
 
 ```
-#option=:other
-allones=1.0 .+0*G.hFacC
-vol=[b(allones) for b in M.h_sum]
+#option=:loops
+M.tmp2d.=M.v_int[1](allones)
+vol=[b(M.tmp2d) for b in M.h_sum]
 ```
 
 or 
 
 ```
-#option=:hv
-M.tmp2d.=M.v_int[1](allones)
-vol=[b(M.tmp2d) for b in M.h_sum]
+#option=:streamlined_loop
+allones=1.0 .+0*G.hFacC
+vol=[b(allones) for b in M.h_sum]
 ```
 """
-function loop_3d(boxes::gridmask; files=String[], var=:THETA, rd=read)
+function loops(mask::gridmask; files=String[], var=:THETA, rd=read)
   nt=length(files)
-  nb=length(boxes.names)
-  BA=SharedArray{Float64}(nb,nt)
+  nh=length(mask.names)
+  nv=length(mask.depths)
+  BA=SharedArray{Float64}(nh,nv,nt)
   @sync @distributed for t in 1:nt
     mod(t,10)==0 ? println(t) : nothing
     F=files[t]
     ext=split(F,".")[end]
-    boxes.tmp3d.=rd(F,var,boxes.tmp3d)
-    BA[:,t]=[b(boxes.tmp3d) for b in boxes.h_sum]
+    mask.tmp3d.=rd(F,var,t,mask.tmp3d)
+    for layer in 1:nv
+       mask.tmp2d.=mask.v_int[layer](mask.tmp3d)
+       BA[:,layer,t]=[b(mask.tmp2d) for b in mask.h_sum]
+    end
     GC.gc()
   end
   BA
 end
 
-##
+"""
+    streamlined_loop(mask::gridmask; files=String[], var=:THETA, rd=read)
 
-function loop_hv(boxes::gridmask; files=String[], var=:THETA, rd=read)
+Alternate approach to loops, where loops are streamlined in a single dimension.
+"""
+function streamlined_loop(mask::gridmask; files=String[], var=:THETA, rd=read)
   nt=length(files)
-  nh=length(boxes.names)
-  nv=length(boxes.depths)
-  BA=SharedArray{Float64}(nh,nv,nt)
+  nb=length(mask.names)
+  BA=SharedArray{Float64}(nb,nt)
   @sync @distributed for t in 1:nt
     mod(t,10)==0 ? println(t) : nothing
     F=files[t]
     ext=split(F,".")[end]
-    boxes.tmp3d.=rd(F,var,boxes.tmp3d)
-    for layer in 1:nv
-       boxes.tmp2d.=boxes.v_int[layer](boxes.tmp3d)
-       BA[:,layer,t]=[b(boxes.tmp2d) for b in boxes.h_sum]
-       #dH0[i]=nansum(write(dT*(G.b_i.==i)))
-    end
+    mask.tmp3d.=rd(F,var,mask.tmp3d)
+    BA[:,t]=[b(mask.tmp3d) for b in mask.h_sum]
     GC.gc()
   end
   BA
 end
 
+##
+
 end
 

src/Interpolation.jl

@@ -247,7 +247,7 @@ function interpolation_setup(;Γ=NamedTuple(),
         lat=[j for i=-179.:2.0:179., j=-89.:2.0:89.])
 
         if isempty(Γ)
-                fil=joinpath(MeshArrays.MA_datadep("interp_halfdeg"),"interp_coeffs_halfdeg.jld2")
+                fil=joinpath(MeshArrays.mydatadep("interp_halfdeg"),"interp_coeffs_halfdeg.jld2")
         else
 		(f,i,j,w)=InterpolationFactors(Γ,vec(lon),vec(lat))
                 fil=tempname()*"_interp_coeffs.jld2"

src/MeshArrays.jl

@@ -50,7 +50,9 @@ export StereographicProjection, isosurface
 #export location_is_out, NeighborTileIndices_dpdo, NeighborTileIndices_cs, RelocationFunctions_cs
 #export update_location_cs!, update_location_llc!, update_location_dpdo!
 
-function MA_datadep end
+function mydatadep end
+
+MA_datadep=mydatadep
 export MA_datadep
 
 function plot_examples end; export plot_examples

src/Types.jl

@@ -256,6 +256,11 @@ end
     gridmask
 
 gridmask data structure.
+
+```
+G,M,files=Integration.example()
+M
+```
 """
 Base.@kwdef struct gridmask
   map::MeshArray

src/demo.jl

@@ -110,21 +110,62 @@ module demo
     """
         ocean_basins()
 
-    Map of ocean basins on ECCO grid.
+    Mask of ocean basins on ECCO grid (LLC90).
 
     ```
     basins=demo.ocean_basins()
+    AtlExt=extended_basin(basins,:Atl)
     ```
     """
     function ocean_basins()        
-        γ=GridSpec("LatLonCap",GRID_LLC90)
+        γ=GridSpec(ID=:LLC90)
         fil_basin=joinpath(GRID_LLC90,"v4_basin.bin")
-        basins=(map=read(fil_basin,MeshArray(γ,Float32)),
+        basins=(mask=read(fil_basin,MeshArray(γ,Float32)),
             name=["Pacific","Atlantic","indian","Arctic","Bering Sea","South China Sea","Gulf of Mexico",
             "Okhotsk Sea","Hudson Bay","Mediterranean Sea","Java Sea","North Sea","Japan Sea",
             "Timor Sea","East China Sea","Red Sea","Gulf","Baffin Bay","GIN Seas","Barents Sea"])
     end
 
+
+    list_Atl=["Atlantic","Gulf of Mexico","Hudson Bay","Mediterranean Sea","North Sea","Baffin Bay","GIN Seas"]
+    list_Pac=["Pacific","Bering Sea","Okhotsk Sea","Japan Sea","East China Sea"]
+    list_Ind=["indian","South China Sea","Java Sea","Timor Sea","Red Sea","Gulf"]
+    list_Arc=["Arctic","Barents Sea"]
+
+    """
+        extended_basin(basins,ID=:Atl)
+
+    Consolidate basins mask to include marginal seas. 
+    
+    note : has only be tested on the ECCO grid (LLC90). 
+
+    ```
+    basins=demo.ocean_basins()
+    AtlExt=extended_basin(basins,:Atl)
+    ```
+    """
+    function extended_basin(basins,ID=:Atl)
+        list_basins=if ID==:Pac
+            list_Pac
+        elseif ID==:Atl
+            list_Atl
+        elseif ID==:Ind
+            list_Ind
+        elseif ID==:Arc
+            list_Arc
+        else
+            error("unknown basin")
+        end
+
+        mask=0*basins.mask
+        for i in list_basins
+            println(i)
+            jj=findall(basins.name.==i)[1]
+            mask.=mask+1.0*(basins.mask.==jj)
+        end
+        mask
+    end
+
     """
         download_polygons(ID::String)
 
@@ -135,9 +176,9 @@ module demo
     """
     function download_polygons(ID::String)
         if ID=="ne_110m_admin_0_countries.shp"
-            joinpath(MeshArrays.MA_datadep("countries_shp1"),"ne_110m_admin_0_countries.shp")
+            joinpath(MeshArrays.mydatadep("countries_shp1"),"ne_110m_admin_0_countries.shp")
         elseif ID=="countries.geojson"
-            joinpath(MeshArrays.MA_datadep("countries_geojson1"),"countries.geojson")
+            joinpath(MeshArrays.mydatadep("countries_geojson1"),"countries.geojson")
         else
             error("unknown data dependency")
         end
@@ -153,7 +194,7 @@ module demo
         lat=[j for i=-0.05:dx:359.95, j=-89.95:dx:89.95]; 
         lon=[i for i=-0.05:dx:359.95, j=-89.95:dx:89.95];
 
-        earth_jpg=joinpath(MeshArrays.MA_datadep("basemap_jpg1"),
+        earth_jpg=joinpath(MeshArrays.mydatadep("basemap_jpg1"),
         "Blue_Marble_Next_Generation_%2B_topography_%2B_bathymetry.jpg")
 
         earth_img=read_JLD2(earth_jpg)

test/runtests.jl

@@ -50,7 +50,7 @@ end
 end
 
 @testset "Regional Integration:" begin
-    G,M,files=Integration.example(option=:hv)
+    G,M,files=Integration.example()
 
     allones=1.0 .+0*G.hFacC
     vol0=sum(G.RAC*G.DRF*G.hFacC)
@@ -59,7 +59,7 @@ end
     vol=[b(M.tmp2d) for b in M.h_sum]
     @test isapprox(sum(vol),vol0)
 
-    G,M,files=Integration.example(option=:full)
+    G,M,files=Integration.example(option=:streamlined_loop)
     vol=[b(allones) for b in M.h_sum]
     @test isapprox(sum(vol),vol0)
 
@@ -230,6 +230,7 @@ end
     λ=interpolation_setup()
 
 	basins=demo.ocean_basins()
+    AtlExt=demo.extended_basin(basins,:Atl)
 	sections,path_sec=demo.ocean_sections(Γ)
     my_section=demo.one_section(Γ,[127 127],[-25 -68])