diff --git a/Project.toml b/Project.toml index 6cd266e..cdbd7fc 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MeshArrays" uuid = "cb8c808f-1acf-59a3-9d2b-6e38d009f683" authors = ["gaelforget "] -version = "0.3.14" +version = "0.3.15" [deps] CatViews = "81a5f4ea-a946-549a-aa7e-2a7f63a27d31" diff --git a/docs/src/API.md b/docs/src/API.md index c90a95a..34feef5 100644 --- a/docs/src/API.md +++ b/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 ``` diff --git a/ext/MeshArraysDataDepsExt.jl b/ext/MeshArraysDataDepsExt.jl index 6ff0a88..19d30f7 100644 --- a/ext/MeshArraysDataDepsExt.jl +++ b/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" diff --git a/src/Grids.jl b/src/Grids.jl index 5fd6794..bc14baf 100644 --- a/src/Grids.jl +++ b/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: diff --git a/src/Integration.jl b/src/Integration.jl index 5b537d7..0d78773 100644 --- a/src/Integration.jl +++ b/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 diff --git a/src/Interpolation.jl b/src/Interpolation.jl index 20a42e8..16dfd07 100644 --- a/src/Interpolation.jl +++ b/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" diff --git a/src/MeshArrays.jl b/src/MeshArrays.jl index 34abaae..0eb43d6 100644 --- a/src/MeshArrays.jl +++ b/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 diff --git a/src/Types.jl b/src/Types.jl index f7685d7..677e8c6 100644 --- a/src/Types.jl +++ b/src/Types.jl @@ -256,6 +256,11 @@ end gridmask gridmask data structure. + +``` +G,M,files=Integration.example() +M +``` """ Base.@kwdef struct gridmask map::MeshArray diff --git a/src/demo.jl b/src/demo.jl index 76c3fcb..e0a75ef 100644 --- a/src/demo.jl +++ b/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) diff --git a/test/runtests.jl b/test/runtests.jl index e6f93b0..95a85f3 100644 --- a/test/runtests.jl +++ b/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])