diff --git a/Project.toml b/Project.toml index 09eae63..458f3d5 100644 --- a/Project.toml +++ b/Project.toml @@ -1,12 +1,13 @@ name = "ECCOtour" uuid = "414fb998-b5a0-4a85-b76f-da39d21445ca" authors = ["G Jake Gebbie"] -version = "0.1.18" +version = "0.1.19" [deps] DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" IsopycnalSurfaces = "29a8b8bc-eb2a-4ada-9a47-be49547c4d59" diff --git a/src/ECCOtour.jl b/src/ECCOtour.jl index 8a12cb6..378e704 100644 --- a/src/ECCOtour.jl +++ b/src/ECCOtour.jl @@ -7,8 +7,10 @@ using Statistics, Distributions, FFTW, #using PyPlot import Statistics.mean, Statistics.std, - Base.maximum, Base.minimum, Base.replace!, - IsopycnalSurfaces.vars2sigma1, IsopycnalSurfaces.sigma1grid + Base.maximum, Base.minimum, Base.replace!, + Base.write, + IsopycnalSurfaces.vars2sigma1, + IsopycnalSurfaces.sigma1grid export hanncoeffs, hannsum, hannsum!, hannfilter export get_filtermatrix, matrixfilter, matrixspray, columnscale! @@ -25,12 +27,9 @@ export regularlatgrid, LLCcropC, LLCcropG, croplimitsLLC export latgridAntarctic, latgridArctic, latgridRegular export timestamp_monthly_v4r4, netcdf2dict, write_vars export mdsio2dict, mdsio2sigma1, ncwritefromtemplate -export netcdf2sigma1, mdsio2regularpoles -export writeregularpoles, vars2regularpoles, var2regularpoles - +export netcdf2sigma1, writeregularpoles export sigma2grid, vars2sigma, mdsio2sigma, mdsio2sigma2 export θbudg2sigma, θbudg2sigma2 - export netcdf2regularpoles, factors4regularpoles export regional_mask, apply_regional_mask!, zero2one!, wrapdist export centerlon!, read_netcdf @@ -41,12 +40,27 @@ export velocity2center, rotate_uv, rotate_velocity! export vars2sigma1, sigma1grid export landmask, land2nan! export cons_offset! - - +export RegularpolesParameters, regularpoles +export grid_attributes, times_ecco +export wet_mask, basin_mask include("HannFilter.jl") include("MatrixFilter.jl") include("SeasonalCycle.jl") +include("basins.jl") + +struct RegularpolesParameters{T<:Real,I<:Integer,NT<: NamedTuple} + λC::StepRangeLen + λG::StepRangeLen + ϕC::Vector{T} + ϕG::Vector{T} + nx::I + ny::I + nyarc::I + nyantarc::I + λarc::NT + λantarc::NT +end """ function sigma2grid() @@ -639,7 +653,8 @@ function factors4regularpoles(γ) nyarc = length(ϕCarc) nyantarc = length(ϕCantarc) - return λC,λG,ϕC,ϕG,nx,ny,nyarc,nyantarc,λarc,λantarc + return RegularpolesParameters(λC,λG,ϕC,ϕG,nx,ny,nyarc,nyantarc,λarc,λantarc) + #return λC,λG,ϕC,ϕG,nx,ny,nyarc,nyantarc,λarc,λantarc end """ @@ -897,29 +912,22 @@ function LLCcropC(gcmfield,γ) # There is a regular grid inside the LLC grid. # Subsample/crop just those points. Put them together correctly. ϕ,λ = latlonG(γ) - jarc,jantarc,ϕarc,ϕantarc = croplimitsLLC(ϕ,λ) + jarc,jantarc,_,_ = croplimitsLLC(ϕ,λ) jarc -= 1 # update for C grid, subtract one here - - regfield = gcmfield[1][:,jantarc:jarc] - regfield = vcat(regfield,gcmfield[2][:,jantarc:jarc]) - for i = 4:5 - tmp = reverse(transpose(gcmfield[i]),dims=2) - regfield = vcat(regfield,tmp[:,jantarc:jarc]) - end - - # handle longitudinal wraparound by hand - wrapval = 218 - xwrap = vcat(wrapval+1:size(regfield,1),1:wrapval) - regfield = regfield[xwrap,:] - return regfield + return joinfields(gcmfield,jarc,jantarc) end function LLCcropG(gcmfield,γ) # There is a regular grid inside the LLC grid. # Subsample/crop just those points. Put them together correctly. ϕ,λ = latlonG(γ) - jarc,jantarc,ϕarc,ϕantarc = croplimitsLLC(ϕ,λ) + jarc,jantarc,_,_ = croplimitsLLC(ϕ,λ) + return joinfields(gcmfield,jarc,jantarc) +end + +function joinfields(gcmfield,jarc,jantarc) + # no tile 3 (i.e., Arctic tile) because it is not regular regfield = gcmfield[1][:,jantarc:jarc] regfield = vcat(regfield,gcmfield[2][:,jantarc:jarc]) for i = 4:5 @@ -930,8 +938,21 @@ function LLCcropG(gcmfield,γ) # handle longitudinal wraparound by hand wrapval = 218 xwrap = vcat(wrapval+1:size(regfield,1),1:wrapval) - regfield = regfield[xwrap,:] - return regfield + return regfield[xwrap,:] + #return regfield + + # Slower code + # # Why no tile 3 (i.e., Arctic tile)? Cannot get Arctic vals through this method. + # regfield = vcat(gcmfield[1][:,jantarc:jarc], + # gcmfield[2][:,jantarc:jarc], + # reverse(transpose(gcmfield[4]),dims=2)[:,jantarc:jarc], + # reverse(transpose(gcmfield[5]),dims=2)[:,jantarc:jarc]) + + # # handle longitudinal wraparound by hand + # wrapval = 218 + # xwrap = vcat(wrapval+1:size(regfield,1),1:wrapval) + # regfield = regfield[xwrap,:] + # return regfield end """ @@ -939,15 +960,22 @@ end print year and month given time index assuming using ECCOv4r4 """ -function timestamp_monthly_v4r4(t) - tstart = 1992 + 1/24 - tend = 2018 - tecco = range(tstart,step=1/12,stop=2018) +function timestamp_monthly_v4r4(t::Integer) + tecco = times_ecco() year = Int(floor(tecco[t])) month = ((t-1)%12)+1 println("year ",year," month ",month) return year,month end +function timestamp_monthly_v4r4(datafile::String) + tsteps_permonth = 730.4 + pattern = r"\.(.*?)\." # regex to extract between periods + tt = round(Int64,parse(Int64,match(pattern,datafile)[1])/tsteps_permonth) + year, month = timestamp_monthly_v4r4(tt) + return year,month +end + +times_ecco() = range(1992 + 1/24,step=1/12,stop=2018) notnanorzero(z) = !iszero(z) && !isnan(z) @@ -1184,7 +1212,7 @@ end - `linearinterp::Logical`: optional keyword argument, do linear interpolation?, default = false - `eos::String`: optional key argument for equation of state, default = "JMD95" """ -function mdsio2sigma1(pathin::String,pathout::String,fileroots::Vector{String},γ,pstdz,sig1grid;splorder=3,linearinterp=false,eos="JMD95") +function mdsio2sigma1(pathin::String,pathout::String,fileroots::Vector{String},γ,pstdz,sig1grid;splorder=3,linearinterp=false,eos="JMD95",writefiles=true) # Read All Variables And Puts Them Into "Vars" Dictionary # ideally would be more generic and Float64 would be covered. @@ -1204,7 +1232,9 @@ function mdsio2sigma1(pathin::String,pathout::String,fileroots::Vector{String}, fileprefix = pathout # use first filename to get timestamp filesuffix = "_on_sigma1"*fileroots[1][14:end]*".data" - write_vars(varsσ,fileprefix,filesuffix) + if writefiles + write_vars(varsσ,fileprefix,filesuffix) + end return varsσ end @@ -1296,23 +1326,23 @@ end function netcdf2regularpoles(ncfilename,ncvarname,γ,nx,ny,nyarc,λarc,nyantarc,λantarc) vars = netcdf2dict(ncfilename,ncvarname,γ) - varsregpoles = vars2regularpoles(vars,γ,nx,ny,nyarc,λarc,nyantarc,λantarc) + varsregpoles = regularpoles(vars,γ,nx,ny,nyarc,λarc,nyantarc,λantarc) end """ function mdsio2regularpoles(pathin,filein,γ,nx,ny,nyarc,λarc,nyantarc,λantarc) """ -function mdsio2regularpoles(pathin,filein,γ,nx,ny,nyarc,λarc,nyantarc,λantarc) - +function regularpoles(pathin,filein,γ,params::RegularpolesParameters) #nx,ny,nyarc,λarc,nyantarc,λantarc) vars = mdsio2dict(pathin,filein,γ) Γ = GridLoad(γ;option="full") rotate_velocity!(vars,Γ) - varsregpoles = vars2regularpoles(vars,γ,nx,ny,nyarc,λarc,nyantarc,λantarc) - + # return vars2regularpoles(vars,γ,params) #nx,ny,nyarc,λarc,nyantarc,λantarc) + return regularpoles(vars,γ,params) #nx,ny,nyarc,λarc,nyantarc,λantarc) end -function vars2regularpoles(vars::Dict{String,MeshArrays.gcmarray{T,2,Matrix{T}}},γ,nx,ny,nyarc,λarc,nyantarc,λantarc) where T<:AbstractFloat +#function vars2regularpoles(vars::Dict{String,MeshArrays.gcmarray{T,2,Matrix{T}}},γ,nx,ny,nyarc,λarc,nyantarc,λantarc) where T<:AbstractFloat +function regularpoles(vars::Dict{String,MeshArrays.gcmarray{T,2,Matrix{T}}},γ,params) where T<:AbstractFloat varsregpoles = Dict{String,Array{T,3}}() NaNT = zero(T)/zero(T) @@ -1324,27 +1354,28 @@ function vars2regularpoles(vars::Dict{String,MeshArrays.gcmarray{T,2,Matrix{T}}} nz = size(varvals,2) #pre-allocate dict - varsregpoles[varname] = fill(NaNT,(nx,ny,nz)) + varsregpoles[varname] = fill(NaNT,(params.nx,params.ny,nz)) for zz = 1:nz # get regular grid by cropping - varsregpoles[varname][:,:,zz]=var2regularpoles(varvals[:,zz],γ,nx,ny,nyarc,λarc,nyantarc,λantarc) + #varsregpoles[varname][:,:,zz]=var2regularpoles(varvals[:,zz],γ,nx,ny,nyarc,λarc,nyantarc,λantarc) + varsregpoles[varname][:,:,zz]=regularpoles(varvals[:,zz],γ,params) #nx,ny,nyarc,λarc,nyantarc,λantarc) end end return varsregpoles end """ -function vars2regularpoles(vars::Dict{String,MeshArrays.gcmarray{T,1,Matrix{T}}},γ,nx,ny,nyarc,λarc,nyantarc,λantarc) +function regularpoles(vars::Dict{String,MeshArrays.gcmarray{T,1,Matrix{T}}},γ,nx,ny,nyarc,λarc,nyantarc,λantarc) variables interpolated onto regularpoles grid """ -function vars2regularpoles(vars::Dict{String,MeshArrays.gcmarray{T,1,Matrix{T}}},γ,nx,ny,nyarc,λarc,nyantarc,λantarc) where T<:AbstractFloat +#function vars2regularpoles(vars::Dict{String,MeshArrays.gcmarray{T,1,Matrix{T}}},γ,nx,ny,nyarc,λarc,nyantarc,λantarc) where T<:AbstractFloat +function regularpoles(vars::Dict{String,MeshArrays.gcmarray{T,1,Matrix{T}}},γ,params) where T<:AbstractFloat varsregpoles = Dict{String,Matrix{T}}() for (varname, varvals) in vars - - varsregpoles[varname] = var2regularpoles(varvals,γ,nx,ny,nyarc,λarc,nyantarc,λantarc) - + #varsregpoles[varname] = var2regularpoles(varvals,γ,nx,ny,nyarc,λarc,nyantarc,λantarc) + varsregpoles[varname] = regularpoles(varvals,γ,params) end return varsregpoles end @@ -1353,7 +1384,8 @@ end var2regularpoles Take one gcmarray in memory, put on regularpoles grid """ -function var2regularpoles(var,γ,nx,ny,nyarc,λarc,nyantarc,λantarc) +#function var2regularpoles(var,γ,nx,ny,nyarc,λarc,nyantarc,λantarc) +function regularpoles(var,γ,params) T = eltype(var) NaNT = zero(T)/zero(T) @@ -1361,28 +1393,33 @@ function var2regularpoles(var,γ,nx,ny,nyarc,λarc,nyantarc,λantarc) # remove contamination from land # this is a problem for masks # going to make the code slow with deepcopy but don't want mutation - varnative = deepcopy(var) - land2nan!(varnative,γ) + varnative = MeshArrays.mask(var,NaN,0.0) + + # older options + #varnative = deepcopy(var) + #varnative = var + #land2nan!(varnative,γ) #pre-allocate output - varregpoles = fill(NaNT,(nx,ny)) + varregpoles = fill(NaNT,(params.nx,params.ny)) # get regular grid by cropping θcrop = LLCcropC(varnative,γ) # interpolate to "regularpoles" - θarc = reginterp(varnative,nx,nyarc,λarc) - θantarc = reginterp(varnative,nx,nyantarc,λantarc) - varregpoles=hcat(θantarc',θcrop,θarc') + θarc = reginterp(varnative,params.nx,params.nyarc,params.λarc) + θantarc = reginterp(varnative,params.nx,params.nyantarc,params.λantarc) + varregpoles = hcat(θantarc',θcrop,θarc') return varregpoles - end """ function land2nan!(msk,γ) - Replace surface land points with NaN +Replace surface land points with NaN + +Deprecate this function: slower than MeshArrays.mask """ function land2nan!(msk,γ) land = landmask(γ) @@ -1403,9 +1440,24 @@ function landmask(γ;level=1) end """ -function writeregularpoles(vars,γ,pathout,filesuffix,filelog,λC,lonatts,ϕC,latatts,z,depthatts) +`function write(vars::Dict{String,Array{Float32,3}}, + params::RegularpolesParameters, + γ::MeshArrays.gcmgrid, + pathout, + filesuffix, + filelog, + gridatts)` + + """ -function writeregularpoles(vars::Dict{String,Array{Float32,3}},γ,pathout,filesuffix,filelog,λC,lonatts,ϕC,latatts,z,depthatts) +function write(vars::Dict{String,Array{Float32,3}}, + params::RegularpolesParameters, + γ::MeshArrays.gcmgrid, + pathout, + filesuffix, + filelog, + gridatts) + #,ϕC,latatts,z,depthatts) for (varname,varvals) in vars println(varname) @@ -1429,17 +1481,18 @@ function writeregularpoles(vars::Dict{String,Array{Float32,3}},γ,pathout,filesu field = varname end - if varname == "p" - fieldDict = Dict("fldname" => "p","title" => "standard pressure", "units" => "dbar", "levs" => length(z)) + if gridatts.depth["longname"] == "Sigma-1" + depthname = "sigma-1" + z = sigma1grid("mixed layer") else - fieldDict = read_available_diagnostics(field,filename=filelog) - + depthname = "depth" + z = depthlevels(γ) end - if depthatts["longname"] == "Sigma-1" - depthname = "sigma-1" + if varname == "p" + fieldDict = Dict("fldname" => "p","title" => "standard pressure", "units" => "dbar", "levs" => length(z)) else - depthname = "depth" + fieldDict = read_available_diagnostics(field,filename=filelog) end # make a directory for this output @@ -1458,14 +1511,14 @@ function writeregularpoles(vars::Dict{String,Array{Float32,3}},γ,pathout,filesu fileout, varname, "lon", - λC, - lonatts, + params.λC, + gridatts.lon, "lat", - ϕC, - latatts, + params.ϕC, + gridatts.lat, depthname, z, - depthatts, + gridatts.depth, atts = varatts, t = NC_FLOAT # Float32 ) @@ -2238,4 +2291,24 @@ function exch_UV_cs3D(fldU::MeshArrays.gcmarray{T, 2, Matrix{T}}, end +grid_attributes() = ( + lon = Dict("longname" => "Longitude", "units" => "degrees east"), + lat = Dict("longname" => "Latitude", "units" => "degrees north"), + depth = Dict("longname" => "Depth", "units" => "m") +) + +# depth is really the vertical coordinate +sigmagrid_attributes() = ( + lon = Dict("longname" => "Longitude", "units" => "degrees east"), + lat = Dict("longname" => "Latitude", "units" => "degrees north"), + depth = Dict("longname" => "Sigma-1", "units" => "kg/m^3 - 1000") +) + +function wet_mask(Γ) + μ =Γ.hFacC[:,1] + μ[findall(μ.>0.0)].=1.0 + μ[findall(μ.==0.0)].=NaN + return μ +end + end #module diff --git a/src/basins.jl b/src/basins.jl new file mode 100644 index 0000000..44609fc --- /dev/null +++ b/src/basins.jl @@ -0,0 +1,148 @@ +basinlist()=["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"] + +""" + function basin_mask(basin_name,γ;hemisphere=nothing,southlat=nothing,northlat=nothing,Lsmooth=nothing) + + Make a mask based on Gael Forget's definitions of ocean basins and sub-basins. + Note: This mask contains float values. It could be more efficient with a BitArray. + +# Arguments +- `basin_name::Union{String,Vector{String}}`: string options are Arctic, Atlantic, Baffin Bay, Barents Sea, Bering Sea, +East China Sea, GIN Seas, Gulf, Gulf of Mexico, Hudson Bay, indian, Japan Sea, Java Sea, +Mediterranean Sea, North Sea, Okhotsk Sea, Pacific, Red Sea, South China Sea, Timor Sea. +-`hemisphere::Symbol`: optional argument with values `:north`, `:south`, and `:both` +-'southlat::Number': optional argument specifying southerly latitude of mask +-'northlat::Number': optional argument specifying northerly latitude of mask +-`Lsmooth::Number`: smoothing lengthscale in grid points +-'southlat::Number': optional argument specifying southerly latitude of mask +-'northlat::Number': optional argument specifying northerly latitude of mask +# Output +- 'mask': space and time field of surface forcing, value of zero inside +designated lat/lon rectangle and fading to 1 outside sponge zone on each edge. This is +because this field ends up being SUBTRACTED from the total forcing +""" +function basin_mask(basin_name) + # open MATLAB files created elsewhere + file = matopen(datadir("basin_grids/GRID_LLC90_"*basin_name)) + for ii in 1:5 + mask[ii] = read(file,basin_name*"_mask"*string(ii)) + end + close(file) + return mask +end +function basin_mask(basin_name::String,γ) + # this version takes one string and returns one mask. + pth = MeshArrays.GRID_LLC90 + Γ = GridLoad(γ;option="full") + basins=read(joinpath(pth,"v4_basin.bin"),MeshArray(γ,Float32)) + basin_list=basinlist() + basinID=findall(basin_list.==basin_name)[1] + basinmask=similar(basins) + for ff in 1:5 + basinmask[ff] .= (basins[ff].==basinID) + end + land2nan!(basinmask,γ) + return basinmask +end +function basin_mask(basin_names::Vector,γ;hemisphere=nothing,southlat=nothing,northlat=nothing,Lsmooth=nothing) + # this is the full version, take a vector of basin names, include optional arguments + pth = MeshArrays.GRID_LLC90 + Γ = GridLoad(γ;option="full") + basins=read(joinpath(pth,"v4_basin.bin"),MeshArray(γ,Float32)) + + mask = 0 * basins # needs NaN on land + for (ii,nn) in enumerate(basin_names) + mask += basin_mask(nn,γ) + end + + if !isnothing(hemisphere) + apply_hemisphere_mask!(mask,hemisphere,γ) + end + + if !isnothing(southlat) && !isnothing(northlat) + apply_latitude_mask!(mask,southlat,northlat,γ) + end + + if !isnothing(Lsmooth) + mask = smooth(mask,Lsmooth,γ) + end + + # change NaNs to zeros. + land2zero!(mask,γ) + return mask +end + +""" + function land2zero!(msk,γ) + +# see tests for different masking from MeshArrays +that could replace this function +""" +function land2zero!(msk,γ) + land = landmask(γ) + for ff in eachindex(msk) + msk[ff][land[ff]] .= zero(eltype(msk)) + end +end + +""" + function smooth(msk::MeshArrays.gcmarray,lengthscale) + + Smooth a gcmarray with a lengthscale of `X` points + + Based off Gael Forget, MeshArrays.jl +""" +function smooth(msk::MeshArrays.gcmarray,X,γ) + Γ = GridLoad(γ;option="full") + DXCsm=X*Γ.DXC; DYCsm=X*Γ.DYC; + #apply smoother + land2nan!(msk,γ) + return msk_smooth=MeshArrays.smooth(msk,DXCsm,DYCsm,Γ); +end + +""" + function apply_hemisphere_mask!(mask,hemisphere,γ) + + overlay a hemispheric mask on a previous `mask` + in-place function + hemisphere options are `:north`,`:south`, and `:both` + + Note: both has not been tested with this version +""" +function apply_hemisphere_mask!(mask,hemisphere,γ) + Γ = GridLoad(γ;option="full") + if hemisphere == :north + hemisphere_mask = Γ.YC .> 0.0; + elseif hemisphere == :south #South + hemisphere_mask = Γ.YC .< 0.0; + elseif hemisphere == :both #both + hemisphere_mask = Γ.YC .> 0.0 || Γ.YC .≤ 0.0; #optional argument? + else + error("no definition for hemisphere") + end + + # need loop to assign/mutate mask + for ff in 1:5 + mask[ff] .*= hemisphere_mask[ff] + end +end + +function apply_latitude_mask!(mask,southlat,northlat,γ) + Γ = GridLoad(γ;option="full") + if southlat < northlat + southcondition = southlat .≤ Γ.YC; + for ff in 1:5 + mask[ff] .*= southcondition[ff] + end + northcondition = Γ.YC .≤ northlat; + for ff in 1:5 + mask[ff] .*= northcondition[ff] + end + else + error("choose a southerly latitude with a value less than the northerly latitude") + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 1f2ebc3..d45a333 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,12 +6,13 @@ using MeshArrays using Statistics using Dierckx using NetCDF +using Downloads #using GoogleDrive @testset "ECCOtour.jl" begin - + pth = MeshArrays.GRID_LLC90 - γ = GridSpec("LatLonCap",pth) + global γ = GridSpec("LatLonCap",pth) Γ = GridLoad(γ;option="full") nf = length(γ.fSize) z = depthlevels(γ) @@ -28,6 +29,10 @@ using NetCDF !ispath(datadir()) && mkdir(datadir()) + @testset "basin mask" begin + include(testdir("test_basinmask.jl")) + end + cd(srcdir()) # workaround: use a shell script #run(`sh $srcdir/download_google_drive.sh`) @@ -35,7 +40,7 @@ using NetCDF # transport > 40 MB, hits virus scanner, move to WHOI website url = "https://www2.whoi.edu/staff/ggebbie/wp-content/uploads/sites/146/2024/03/trsp_3d_set1.0000000732.tar_.gz" - download(url,"trsp_3d_set1.0000000732.tar.gz") + Downloads.download(url,"trsp_3d_set1.0000000732.tar.gz") run(`tar xvzf state_3d_set1.0000000732.tar.gz`) run(`tar xvzf trsp_3d_set1.0000000732.tar.gz`) @@ -114,26 +119,28 @@ using NetCDF # Set up Cartesian grid for interpolation. # Time for a structure. - λC,λG,ϕC,ϕG,nx,ny,nyarc,nyantarc,λarc,λantarc = - factors4regularpoles(γ) + rp_params = factors4regularpoles(γ) @testset "regularpoles 3d state" begin filein = fileroots[1] pathin = datadir() - @time varsregpoles = mdsio2regularpoles(pathin,filein,γ,nx,ny,nyarc,λarc,nyantarc,λantarc) + @time varsregpoles = regularpoles(pathin,filein,γ,rp_params) #nx,ny,nyarc,λarc,nyantarc,λantarc) filesuffix = "suffix.nc" pathout = pathin filelog = srcdir("available_diagnostics.log") - lonatts = Dict("longname" => "Longitude", "units" => "degrees east") - latatts = Dict("longname" => "Latitude", "units" => "degrees north") - depthatts = Dict("longname" => "Depth", "units" => "m") + gridatts = grid_attributes() - @time writeregularpoles(varsregpoles,γ,pathout,filesuffix,filelog,λC,lonatts, - ϕC,latatts,z,depthatts) - + @time write(varsregpoles, + rp_params, + γ, + pathout, + filesuffix, + filelog, + gridatts) + @test maximum(filter(!isnan,varsregpoles["SALT"])) < 50. @test minimum(filter(!isnan,varsregpoles["SALT"])) > 0. @@ -147,7 +154,7 @@ using NetCDF filein = fileroots[2] pathin = datadir() - @time varsregpoles = mdsio2regularpoles(pathin,filein,γ,nx,ny,nyarc,λarc,nyantarc,λantarc) + @time varsregpoles = regularpoles(pathin,filein,γ,rp_params) @test maximum(filter(!isnan,varsregpoles["NVELMASS"])) < 2. @test minimum(filter(!isnan,varsregpoles["NVELMASS"])) > -2. @@ -156,13 +163,11 @@ using NetCDF @testset "regularpoles dxc" begin # load centered longitude - # dxc = γ.read(γ.path*"XC.data",MeshArray(γ,Float64)) vars = Dict("XC" => γ.read(γ.path*"XC.data",MeshArray(γ,Float64))) - dxc_regpoles = vars2regularpoles(vars,γ,nx,ny,nyarc,λarc,nyantarc,λantarc) - + @time dxc_regpoles = ECCOtour.regularpoles(vars,γ,rp_params) yy = 100 for xx = 1:50 - @test isapprox(dxc_regpoles["XC"][xx,yy],λC[xx], rtol=1e-6) + @test isapprox(dxc_regpoles["XC"][xx,yy],rp_params.λC[xx], rtol=1e-6) end end end @@ -202,7 +207,6 @@ using NetCDF end - @testset "MeshArrays" begin ###################################### # Test the statistics for MeshArrays. diff --git a/test/test_basinmask.jl b/test/test_basinmask.jl new file mode 100644 index 0000000..a937390 --- /dev/null +++ b/test/test_basinmask.jl @@ -0,0 +1,17 @@ +maskname = ["Pacific","South China Sea","East China Sea","Okhotsk Sea","Java Sea","Japan Sea","Timor Sea"] +hemisphere = :both +Lsmooth = 5 +southlat = -15 +northlat = 15 + +msk = basin_mask(maskname,γ,southlat=southlat,northlat=northlat) +@test maximum(msk) == 1.0 +@test minimum(msk) == 0.0 + +@test iszero(maximum(nancount(msk))) +@test sum(nancount(MeshArrays.mask(msk,NaN,0.0))) > 0 + +msk_smooth = basin_mask(maskname,γ,southlat=southlat,northlat=northlat,Lsmooth=Lsmooth) +maximum(nancount(msk_smooth)) +# any values between 0 and 1? yes. +@test maximum(MeshArrays.mask(msk,0.0)) > maximum(MeshArrays.mask(msk_smooth,0.0)) # smooth always less than 1