diff --git a/NEWS.md b/NEWS.md new file mode 100644 index 00000000..f0a3d25d --- /dev/null +++ b/NEWS.md @@ -0,0 +1,11 @@ +# News + +## v0.2 + +* a list of pre-defined lon-lat boxes was added that can be accessed when calling `getCubeData`, so for example `getCubeData(c,region="Germany")` will extract a Germany from the Cube. +* Support for masks was extended, this includes multiple changes: + - every cube now has a `properties` field, which contains a `Dict{String,Any}` where additional cube propertyies can be stored + - if an Integer-valued cube is used as a categorical mask one can add a Dict containg the mapping from value to label name to the cubes properties, eg `c.properties["labels"]=Dict(1=>"low",2=>"high")` + - `plotMAP` will respect the defined label properties and use the labels for its legend + - if a labeled cube is used in the `by` argument when fitting an `OnlineStat`, the output axis is automatically created + - accessing a static variable (currently only water_mask or country_mask) through `getCubeData` will by default only return the first time step diff --git a/docs/src/analysis.md b/docs/src/analysis.md index a38b0015..13c90477 100644 --- a/docs/src/analysis.md +++ b/docs/src/analysis.md @@ -82,10 +82,11 @@ statistics by a certain axis, it has to be specified using the `by` keyword argu `by` accepts a vector of axes types and up to one datacube that can serve as a mask. If such a data cube is supplied, the statistics are split by the unique values in the mask. One can pass a function `cfun` that transforms the mask values into an index in the range `1..N` that defines the -index where the new value is going to be put to. If a mask is supplied, one also needs to supply an +index where the new value is going to be put to. If a mask is supplied, it must have either a `labels` property, +which is a `Dict{T,String}` mapping the numerical mask value to the value name. Alternatively on can supply an `outAxis` argument that describes the resulting output axis. -This all gets clearer with a little example. suppose we want to calculate the mean of GPP, NEE and TER +This all gets clearer with two small examples. suppose we want to calculate the mean of GPP, NEE and TER under the condition that Tair<280K and Tair>280K over all time steps and grid cells. This is achieved through the following lines of code: @@ -129,3 +130,17 @@ b=IOBuffer() show(b,MIME"text/html"(),p) Documenter.Documents.RawHTML(takebuf_string(b)) ``` + +A second example would be that we want to calculate averages of the fluxes according to +a country mask. + +```julia +import OnlineStats +vars = ["gross_primary_productivity","net_ecosystem_exchange","terrestrial_ecosystem_respiration"] +m = getCubeData(ds,variable="country_mask",longitude=lons,latitude=lats) +cube = getCubeData(ds,variable=vars,longitude=lons,latitude=lats) + +mT = mapCube(OnlineStats.Mean,cube,by=[m,VariableAxis], cfun=splitTemp, outAxis=outAxis) +``` + +This will split the cube by country and variable and compute averages over the input variables. diff --git a/src/CubeAPI/CachedArrays.jl b/src/CubeAPI/CachedArrays.jl index 12ab33c4..27b8719f 100644 --- a/src/CubeAPI/CachedArrays.jl +++ b/src/CubeAPI/CachedArrays.jl @@ -337,12 +337,18 @@ write_subblock!{T,N}(x::MaskedCacheBlock{T,N},y::Any,block_size::CartesianIndex{ function write_subblock!{T,N}(x::MaskedCacheBlock{T,N},y::TempCube{T,N},block_size::CartesianIndex{N}) filename=joinpath(y.folder,tofilename(x.position)) - #println("Writing to file $filename") ncwrite(x.data,filename,"cube") ncwrite(x.mask,filename,"mask") ncclose(filename) x.iswritten=false end +function write_subblock!{T}(x::MaskedCacheBlock{T,0},y::TempCube{T,0},block_size::CartesianIndex{0}) + filename=joinpath(y.folder,tofilename(x.position)) + ncwrite([x.data[1]],filename,"cube") + ncwrite([x.mask[1]],filename,"mask") + ncclose(filename) + x.iswritten=false +end function read_subblock!{T,N}(x::MaskedCacheBlock{T,N},y::TempCube{T,N},block_size::CartesianIndex{N}) if block_size==y.block_size filename=joinpath(y.folder,tofilename(x.position)) @@ -356,6 +362,19 @@ function read_subblock!{T,N}(x::MaskedCacheBlock{T,N},y::TempCube{T,N},block_siz end end +function read_subblock!{T}(x::MaskedCacheBlock{T,0},y::TempCube{T,0},block_size::CartesianIndex{0}) + filename=joinpath(y.folder,tofilename(x.position)) + vc=NetCDF.open(filename,"cube") + vm=NetCDF.open(filename,"mask") + x.data[1]=vc[1] + x.mask[1]=vm[1] + ncclose(filename) +end +function getSubRange{T,S<:MaskedCacheBlock}(c::CachedArray{T,0,S};write=false) + c.currentblocks[1].iswritten=write + (c.currentblocks[1].data,c.currentblocks[1].mask) +end + function sync(c::CachedArray) for b in c.currentblocks if b.iswritten diff --git a/src/CubeAPI/CubeAPI.jl b/src/CubeAPI/CubeAPI.jl index 090edd12..18d9e978 100644 --- a/src/CubeAPI/CubeAPI.jl +++ b/src/CubeAPI/CubeAPI.jl @@ -13,6 +13,7 @@ importall .Mask using DataStructures using Base.Dates import Requests.get +import DataStructures.OrderedDict using LightXML type ConfigEntry{LHS} @@ -20,7 +21,10 @@ type ConfigEntry{LHS} rhs end - +#This is probably not the final solution, we define a set of static variables that are treated differently when reading +const static_vars = Set(["water_mask","country_mask"]) +const countrylabels = include("countrylabels.jl") +const known_labels = Dict("water_mask"=>Dict(0x01=>"land",0x02=>"water"),"country_mask"=>countrylabels) " A data cube's static configuration information. @@ -228,6 +232,7 @@ immutable SubCube{T,C} <: AbstractSubCube{T,3} lonAxis::LonAxis latAxis::LatAxis timeAxis::TimeAxis + properties::Dict{String} end axes(s::SubCube)=CubeAxis[s.lonAxis,s.latAxis,s.timeAxis] @@ -287,6 +292,7 @@ immutable SubCubeV{T,C} <: AbstractSubCube{T,4} latAxis::LatAxis timeAxis::TimeAxis varAxis::VariableAxis + properties::Dict{String} end """ @@ -313,6 +319,54 @@ s.perm[4]==1 ? length(s.parent.lonAxis) : s.perm[4]==2 ? length(s.parent.latAxis Base.permutedims{T}(c::SubCube{T},perm::NTuple{3,Int})=SubCubePerm(c,perm) Base.permutedims{T}(c::SubCubeV{T},perm::NTuple{4,Int})=SubCubeVPerm(c,perm) +""" + immutable SubCubeStatic{T, C} <: AbstractCubeData{T,2} + +A view into the data cube with a single static variable. Returned by the `mapCube` +function. + +### Fields + +* `cube::C` Parent cube +* `variable` list of selected variables +* `sub_grid` representation of the subgrid indices +* `time` representation of the selected time steps +* `lonAxis` +* `latAxis` + +""" +immutable SubCubeStatic{T,C} <: AbstractSubCube{T,2} + cube::C #Parent cube + variable::String #Variable + sub_grid::Tuple{Int,Int,Int,Int} #grid_y1,grid_y2,grid_x1,grid_x2 + sub_times::NTuple{6,Int} #y1,i1,y2,i2,ntime,NpY + lonAxis::LonAxis + latAxis::LatAxis + properties::Dict{String} +end + +""" + immutable SubCubeStaticPerm{T} <: AbstractCubeData{T,2} + +Representation of a `SubCubeStatic` with permuted dimensions. +""" +immutable SubCubeStaticPerm{T} <: AbstractSubCube{T,2} + parent::SubCubeV{T} +end +SubCubeStaticPerm{T}(p::SubCubeStatic{T})=SubCubeVPerm{T}(p) +axes(s::SubCubeStatic)=CubeAxis[s.lonAxis,s.latAxis] +axes(s::SubCubeStaticPerm)=CubeAxis[s.parent.latAxis,s.parent.lonAxis] +Base.ndims(s::SubCubeStatic)=2 +Base.ndims(s::SubCubeStaticPerm)=2 +Base.size(s::SubCubeStatic)=(length(s.lonAxis),length(s.latAxis)) +Base.size(s::SubCubeStaticPerm)=(length(s.latAxis),length(s.lonAxis)) + +sgetperm(s::Union{SubCubePerm,SubCubeVPerm})=s.perm +sgetiperm(s::Union{SubCubePerm,SubCubeVPerm})=s.iperm +getperm(s::SubCubeStaticPerm)=(2,1) +getiperm(s::SubCubeStaticPerm)=(2,1) +Base.permutedims{T}(c::SubCubeStatic{T},perm::NTuple{2,Int})=perm==(2,1) ? SubCubeStaticPerm(c) : error("There is only one permutation of a lon-lat cube, so perm must be (2,1)") + """ @@ -333,7 +387,7 @@ Returns a `SubCube` object which represents a view into the original data cube. function getCubeData(cube::UCube;variable=Int[],time=[],latitude=[],longitude=[],region=[]) #First fill empty inputs isempty(variable) && (variable = defaultvariable(cube)) - isempty(time) && (time = defaulttime(cube)) + time==Int[] && (time = defaulttime(cube,variable)) if !isempty(region) haskey(known_regions,region) || error("Region $region not recognized as a known place") ll = known_regions[region] @@ -342,13 +396,15 @@ function getCubeData(cube::UCube;variable=Int[],time=[],latitude=[],longitude=[] end isempty(latitude) && (latitude = defaultlatitude(cube)) isempty(longitude) && (longitude= defaultlongitude(cube)) - getCubeData(cube,variable,time,latitude,longitude) + getCubeData(cube,variable,time,longitude,latitude) end -defaulttime(cube::UCube)=cube.config.start_time,cube.config.end_time-Day(1) -defaultvariable(cube::UCube)=cube.dataset_files +defaulttime(cube::UCube,v)=all_mask(v) ? cube.config.start_time : (cube.config.start_time,cube.config.end_time-Day(1)) +defaultvariable(cube::UCube)=filter(i->i ∉ static_vars,cube.dataset_files) defaultlatitude(cube::UCube)=(-90.0,90.0) defaultlongitude(cube::UCube)=(-180.0,180.0) +all_mask(x::String)=x ∈ static_vars +all_mask(x::Vector{String})=all(all_mask,x) using NetCDF vartype{T,N}(v::NcVar{T,N})=T @@ -372,17 +428,17 @@ function getTimeRanges(c::UCube,y1,y2,i1,i2) end #Convert single input to vectors -function getCubeData{T<:Union{Integer,AbstractString}}(cube::UCube, - variable::Union{AbstractString,Integer,AbstractVector{T}}, - time::Union{Tuple{TimeType,TimeType},TimeType}, - latitude::Union{Tuple{Real,Real},Real}, - longitude::Union{Tuple{Real,Real},Real}) +function getCubeData(cube::UCube, + variable, + time, + longitude, + latitude) isa(time,TimeType) && (time=(time,time)) isa(latitude,Real) && (latitude=(latitude,latitude)) isa(longitude,Real) && (longitude=(longitude,longitude)) - isa(variable,AbstractVector) && isa(eltype(variable),Integer) && (variable=[cube.config.dataset_files[i] for i in variable]) - isa(variable,Integer) && (variable=dataset_files[i]) + !isa(variable,Vector) && (variable=[variable]) + isa(eltype(variable),Integer) && (variable=[cube.config.dataset_files[i] for i in variable]) getCubeData(cube,variable,time,longitude,latitude) end @@ -417,6 +473,7 @@ function getMaskFile(cube::RemoteCube) end end +getLandSeaMask!(mask::Array{UInt8,2},cube::UCube,grid_x1,nx,grid_y1,ny)=getLandSeaMask!(reshape(mask,(size(mask,1),size(mask,2),1)),cube,grid_x1,nx,grid_y1,ny) function getLandSeaMask!(mask::Array{UInt8,3},cube::UCube,grid_x1,nx,grid_y1,ny) filename=getMaskFile(cube) if !isempty(filename) @@ -459,44 +516,21 @@ function getvartype(cube::RemoteCube,variable) eltype(NetCDF.open(datafile,variable)) end - -function getCubeData(cube::UCube, - variable::AbstractString, - time::Tuple{TimeType,TimeType}, - latitude::Tuple{Real,Real}, - longitude::Tuple{Real,Real}) - # This function is doing the actual reading - config=cube.config - - longitude[1]>longitude[2] && throw(ArgumentError("Longitudes must be passed in West-to-East order.")) - - latitude[1]>latitude[2] && (latitude=(latitude[2],latitude[1])) - - grid_y1,grid_y2,grid_x1,grid_x2 = getLonLatsToRead(config,longitude,latitude) - y1,i1,y2,i2,ntime,NpY = getTimesToRead(time[1],time[2],config) - - t=getvartype(cube,variable) - - return SubCube{t,typeof(cube)}(cube,variable, - (grid_y1,grid_y2,grid_x1,grid_x2), - (y1,i1,y2,i2,ntime,NpY), - LonAxis(x2lon(grid_x1,config):config.spatial_res:x2lon(grid_x2,config)), - LatAxis(y2lat(grid_y1,config):-config.spatial_res:y2lat(grid_y2,config)), - TimeAxis(getTimeRanges(cube,y1,y2,i1,i2))) -end - function getCubeData{T<:AbstractString}(cube::UCube, variable::Vector{T}, time::Tuple{TimeType,TimeType}, - latitude::Tuple{Real,Real}, - longitude::Tuple{Real,Real}) + longitude::Tuple{Real,Real}, + latitude::Tuple{Real,Real}) config=cube.config - longitude[1]>longitude[2] && throw(ArgumentError("Longitudes must be passed in West-to-East order.")) + longitude[1]>longitude[2] && throw(ArgumentError("Longitudes $longitude must be passed in West-to-East order.")) + + latitude[1]>latitude[2] && (latitude=(latitude[2],latitude[1])) grid_y1,grid_y2,grid_x1,grid_x2 = getLonLatsToRead(config,longitude,latitude) y1,i1,y2,i2,ntime,NpY = getTimesToRead(time[1],time[2],config) + variableNew=String[] varTypes=DataType[] for i=1:length(variable) @@ -508,14 +542,45 @@ function getCubeData{T<:AbstractString}(cube::UCube, warn("Skipping variable $(variable[i]), not found in Datacube") end end - tnew=reduce(promote_type,varTypes[1],varTypes) - return SubCubeV{tnew,typeof(cube)}(cube,variable, - (grid_y1,grid_y2,grid_x1,grid_x2), - (y1,i1,y2,i2,ntime,NpY), - LonAxis(x2lon(grid_x1,config):config.spatial_res:x2lon(grid_x2,config)), - LatAxis(y2lat(grid_y1,config):-config.spatial_res:y2lat(grid_y2,config)), - TimeAxis(getTimeRanges(cube,y1,y2,i1,i2)), - VariableAxis(variableNew)) + t=reduce(promote_type,varTypes[1],varTypes) + properties=Dict{String,Any}() + if length(variableNew)==1 + if time[1]==time[2] + # This is a static cube and probably small enough to be in memory + c=SubCubeStatic{t,typeof(cube)}(cube,variable[1], + (grid_y1,grid_y2,grid_x1,grid_x2), + (y1,i1,y2,i2,ntime,NpY), + LonAxis(x2lon(grid_x1,config):config.spatial_res:x2lon(grid_x2,config)), + LatAxis(y2lat(grid_y1,config):-config.spatial_res:y2lat(grid_y2,config)), + properties) + d=readCubeData(c) + if haskey(known_labels,variable[1]) + rem_keys = unique(d.data) + all_labels = known_labels[variable[1]] + left_labels = OrderedDict((k,all_labels[k]) for k in rem_keys if k0 - @assert nt<=size(outar,3) - @assert ny<=size(outar,2) - @assert nx<=size(outar,1) + @assert nt==size(outar,3) + @assert ny==size(outar,2) + @assert nx==size(outar,1) @assert size(outar)==size(mask) v=try NetCDF.open(filename,variable); catch - @inbounds for k=1:nt,j=1:ny,i=1:nx - mask[i,j,k]=(mask[i,j,k] | OUTOFPERIOD) - outar[i,j,k]=nanval - end + mask[:] = OUTOFPERIOD + outar[:] = nanval return false end scalefac::T = convert(T,get(v.atts,"scale_factor",one(T))) offset::T = convert(T,get(v.atts,"add_offset",zero(T))) - outar[1:nx,1:ny,1:nt]=v[xr,yr,tstart:(tstart+nt-1)] + NetCDF.readvar!(v,outar,start=[grid_x1,grid_y1,tstart],count=[nx,ny,nt]) missval::T=convert(T,ncgetatt(filename,variable,"_FillValue")) - @inbounds for k=1:nt,j=1:ny,i=1:nx - if (outar[i,j,k] == missval) || isnan(outar[i,j,k]) - mask[i,j,k]=mask[i,j,k] | MISSING - outar[i,j,k]=nanval + @inbounds for i in eachindex(outar) + if (outar[i] == missval) || isnan(outar[i]) + mask[i] = MISSING + outar[i] = nanval else - outar[i,j,k]=outar[i,j,k]*scalefac+offset + outar[i]=outar[i]*scalefac+offset end end ncclose(filename) return true end +gettoffsnt(::AbstractSubCube,r::CartesianRange)=(r.start.I[3] - 1,r.stop.I[3] - r.start.I[3]+1) +gettoffsnt(::SubCubeStatic,r::CartesianRange{CartesianIndex{2}})=(0,1) function _read{T}(s::AbstractSubCube{T},t::NTuple{2},r::CartesianRange) #;xoffs::Int=0,yoffs::Int=0,toffs::Int=0,voffs::Int=0,nx::Int=size(outar,1),ny::Int=size(outar,2),nt::Int=size(outar,3),nv::Int=length(s.variable)) @@ -618,8 +699,7 @@ end nx = r.stop.I[1] - r.start.I[1] + 1 grid_y1 = grid_y1 + r.start.I[2] - 1 ny = r.stop.I[2] - r.start.I[2] +1 - toffs = r.start.I[3] - 1 - nt = r.stop.I[3] - r.start.I[3]+1 + toffs,nt= gettoffsnt(s,r) if toffs > 0 i1 = i1 + toffs if i1 > NpY @@ -652,6 +732,16 @@ end end end + function readAllyears(s::SubCubeStatic,outar,mask,y1,i1,grid_x1,nx,grid_y1,ny,nt,voffs,nv,NpY) + ycur=y1 #Current year to read + i1cur=i1 #Current time step in year + itcur=1 #Current time step in output file + fin = false + while !fin + fin,ycur,i1cur,itcur = readFromDataYear(s.cube,reshape(outar,(size(outar,1),size(outar,2),1)),reshape(mask,(size(mask,1),size(mask,2),1)),s.variable,ycur,grid_x1,nx,grid_y1,ny,itcur,i1cur,nt,NpY) + end + end + function readAllyears(s::SubCubeV,outar,mask,y1,i1,grid_x1,nx,grid_y1,ny,nt,voffs,nv,NpY) for iv in (voffs+1):(nv+voffs) outar2=view(outar,:,:,:,iv-voffs) @@ -684,6 +774,7 @@ function getremFileName(cube::Cube,variable::String) end showVarInfo(cube::SubCube)=showVarInfo(cube,cube.variable) +showVarInfo(cube::SubCubeStatic)=showVarInfo(cube,cube.variable) showVarInfo(cube::SubCubeV)=[showVarInfo(cube,v) for v in cube.variable] function showVarInfo(cube, variable::String) filename=getremFileName(cube.cube,variable) @@ -720,34 +811,37 @@ end show(io::IO,::MIME"text/markdown",v::Vector{CABLABVarInfo})=foreach(x->show(io,MIME"text/markdown"(),x),v) - function readFromDataYear{T}(cube::Cube,outar::AbstractArray{T,3},mask::AbstractArray{UInt8,3},variable,y,grid_x1,nx,grid_y1,ny,itcur,i1cur,ntime,NpY) - filename=joinpath(cube.base_dir,"data",variable,string(y,"_",variable,".nc")) - ntleft = ntime - itcur + 1 - nt = min(NpY-i1cur+1,ntleft) - xr = grid_x1:(grid_x1+nx-1) - yr = grid_y1:(grid_y1+ny-1) - nanval=convert(T,NaN) - #Make some assertions for inbounds - @assert itcur>0 - @assert (itcur+nt-1)<=size(outar,3) - @assert ny<=size(outar,2) - @assert nx<=size(outar,1) - @assert size(outar)==size(mask) - if isfile(filename) - v=NetCDF.open(filename,variable); - scalefac::T = convert(T,get(v.atts,"scale_factor",one(T))) - offset::T = convert(T,get(v.atts,"add_offset",zero(T))) - outar[1:nx,1:ny,itcur:(itcur+nt-1)]=v[xr,yr,i1cur:(i1cur+nt-1)] - missval::T=convert(T,ncgetatt(filename,variable,"_FillValue")) - @inbounds for k=itcur:(itcur+nt-1),j=1:ny,i=1:nx - if (outar[i,j,k] == missval) || isnan(outar[i,j,k]) - mask[i,j,k]=mask[i,j,k] | MISSING - outar[i,j,k]=nanval - else - outar[i,j,k]=outar[i,j,k]*scalefac+offset - end +getNanVal{T<:AbstractFloat}(::Type{T}) = convert(T,NaN) +getNanVal{T<:Integer}(::Type{T}) = typemax(T) + +function readFromDataYear{T}(cube::Cube,outar::AbstractArray{T,3},mask::AbstractArray{UInt8,3},variable,y,grid_x1,nx,grid_y1,ny,itcur,i1cur,ntime,NpY) + filename=joinpath(cube.base_dir,"data",variable,string(y,"_",variable,".nc")) + ntleft = ntime - itcur + 1 + nt = min(NpY-i1cur+1,ntleft) + xr = grid_x1:(grid_x1+nx-1) + yr = grid_y1:(grid_y1+ny-1) + nanval=getNanVal(T) + #Make some assertions for inbounds + @assert itcur>0 + @assert (itcur+nt-1)<=size(outar,3) + @assert ny==size(outar,2) + @assert nx==size(outar,1) + @assert size(outar)==size(mask) + if isfile(filename) + v=NetCDF.open(filename,variable) + scalefac::T = convert(T,get(v.atts,"scale_factor",one(T))) + offset::T = convert(T,get(v.atts,"add_offset",zero(T))) + NetCDF.readvar!(v,view(outar,:,:,itcur:(itcur+nt-1)),start=[grid_x1,grid_y1,i1cur],count=[nx,ny,nt]) + missval::T=convert(T,ncgetatt(filename,variable,"_FillValue")) + @inbounds for k=itcur:(itcur+nt-1),j=1:ny,i=1:nx + if (outar[i,j,k] == missval) || isnan(outar[i,j,k]) + mask[i,j,k]=mask[i,j,k] | MISSING + outar[i,j,k]=nanval + else + outar[i,j,k]=outar[i,j,k]*scalefac+offset end - ncclose(filename) + end + ncclose(filename) else @inbounds for k=itcur:(itcur+nt-1),j=1:ny,i=1:nx mask[i,j,k]=(mask[i,j,k] | OUTOFPERIOD) diff --git a/src/CubeAPI/countrydict.jl b/src/CubeAPI/countrydict.jl index 2c27310d..74b3d785 100644 --- a/src/CubeAPI/countrydict.jl +++ b/src/CubeAPI/countrydict.jl @@ -1,4 +1,15 @@ """ + +### List of Continents + +* Africa +* Asia +* Australia +* Europe +* North America +* South America + + ### List of SREX regions | Short Name | Long Name | @@ -296,6 +307,13 @@ |ZWE|Zimbabwe """ const known_regions = Dict{String,NTuple{4,Float64}}( +"Africa"=>(-17.0,-40.0,51.0,40.0), +"Asia"=>(30.0,5.0,180.0,178.0), +"Australia"=>(110.0,-40.0,155.0,-8.0), +"Europe"=>(-10.0,35.0,33.0,70.0), +"North America"=>(-168.0,10.0,-55.0,80.0), +"South America"=>(-85.0,-60.0,-33.0,15.0), + "Alaska/N.W. Canada" => (-168.0,60.0,-105.0,72.6), "ALA" => (-168.0,60.0,-105.0,72.6), "Amazon" => (-79.7,-20.0,-50.0,11.4), diff --git a/src/CubeAPI/countrylabels.jl b/src/CubeAPI/countrylabels.jl new file mode 100644 index 00000000..b19ba598 --- /dev/null +++ b/src/CubeAPI/countrylabels.jl @@ -0,0 +1,253 @@ +Dict{Int32,String}( +533 => "Aruba", +4 => "Afghanistan", +24 => "Angola", +660 => "Anguilla", +8 => "Albania", +248 => "Aland", +20 => "Andorra", +784 => "United Arab Emirates", +32 => "Argentina", +51 => "Armenia", +16 => "American Samoa", +10 => "Antarctica", +36 => "Ashmore and Cartier Is.", +260 => "Fr. S. and Antarctic Lands", +28 => "Antigua and Barb.", +36 => "Australia", +40 => "Austria", +31 => "Azerbaijan", +108 => "Burundi", +56 => "Belgium", +204 => "Benin", +854 => "Burkina Faso", +50 => "Bangladesh", +100 => "Bulgaria", +48 => "Bahrain", +44 => "Bahamas", +70 => "Bosnia and Herz.", +652 => "St. Barthelemy", +112 => "Belarus", +84 => "Belize", +60 => "Bermuda", +68 => "Bolivia", +76 => "Brazil", +52 => "Barbados", +96 => "Brunei", +64 => "Bhutan", +72 => "Botswana", +140 => "Central African Rep.", +124 => "Canada", +756 => "Switzerland", +152 => "Chile", +156 => "China", +384 => "Ivory Coast", +-99 => "Clipperton I.", +120 => "Cameroon", +180 => "Congo (Kinshasa)", +178 => "Congo (Brazzaville)", +184 => "Cook Is.", +170 => "Colombia", +174 => "Comoros", +132 => "Cape Verde", +188 => "Costa Rica", +-99 => "Coral Sea Is.", +192 => "Cuba", +531 => "Curacao", +136 => "Cayman Is.", +-99 => "N. Cyprus", +196 => "Cyprus", +203 => "Czech Rep.", +276 => "Germany", +262 => "Djibouti", +212 => "Dominica", +208 => "Denmark", +214 => "Dominican Rep.", +12 => "Algeria", +218 => "Ecuador", +818 => "Egypt", +232 => "Eritrea", +-99 => "Dhekelia", +724 => "Spain", +233 => "Estonia", +231 => "Ethiopia", +246 => "Finland", +242 => "Fiji", +238 => "Falkland Is.", +250 => "France", +234 => "Faroe Is.", +583 => "Micronesia", +266 => "Gabon", +275 => "Gaza", +826 => "United Kingdom", +268 => "Georgia", +831 => "Guernsey", +288 => "Ghana", +292 => "Gibraltar", +324 => "Guinea", +270 => "Gambia", +624 => "Guinea Bissau", +226 => "Eq. Guinea", +300 => "Greece", +308 => "Grenada", +304 => "Greenland", +320 => "Guatemala", +316 => "Guam", +328 => "Guyana", +344 => "Hong Kong", +334 => "Heard I. and McDonald Is.", +340 => "Honduras", +191 => "Croatia", +332 => "Haiti", +348 => "Hungary", +360 => "Indonesia", +833 => "Isle of Man", +356 => "India", +-99 => "Indian Ocean Ter.", +86 => "Br. Indian Ocean Ter.", +372 => "Ireland", +364 => "Iran", +368 => "Iraq", +352 => "Iceland", +376 => "Israel", +380 => "Italy", +388 => "Jamaica", +832 => "Jersey", +400 => "Jordan", +392 => "Japan", +-99 => "Baykonur", +398 => "Kazakhstan", +404 => "Kenya", +417 => "Kyrgyzstan", +116 => "Cambodia", +296 => "Kiribati", +659 => "St. Kitts and Nevis", +410 => "S. Korea", +-99 => "Kosovo", +414 => "Kuwait", +418 => "Laos", +422 => "Lebanon", +430 => "Liberia", +434 => "Libya", +662 => "Saint Lucia", +438 => "Liechtenstein", +144 => "Sri Lanka", +426 => "Lesotho", +440 => "Lithuania", +442 => "Luxembourg", +428 => "Latvia", +446 => "Macau", +663 => "St. Martin", +504 => "Morocco", +492 => "Monaco", +498 => "Moldova", +450 => "Madagascar", +462 => "Maldives", +484 => "Mexico", +584 => "Marshall Is.", +807 => "Macedonia", +466 => "Mali", +470 => "Malta", +104 => "Myanmar", +499 => "Montenegro", +496 => "Mongolia", +580 => "N. Mariana Is.", +508 => "Mozambique", +478 => "Mauritania", +500 => "Montserrat", +480 => "Mauritius", +454 => "Malawi", +458 => "Malaysia", +516 => "Namibia", +540 => "New Caledonia", +562 => "Niger", +574 => "Norfolk Island", +566 => "Nigeria", +558 => "Nicaragua", +570 => "Niue", +528 => "Netherlands", +578 => "Norway", +524 => "Nepal", +520 => "Nauru", +554 => "New Zealand", +512 => "Oman", +586 => "Pakistan", +591 => "Panama", +612 => "Pitcairn Is.", +604 => "Peru", +608 => "Philippines", +585 => "Palau", +598 => "Papua New Guinea", +616 => "Poland", +630 => "Puerto Rico", +408 => "N. Korea", +620 => "Portugal", +600 => "Paraguay", +258 => "Fr. Polynesia", +634 => "Qatar", +642 => "Romania", +643 => "Russia", +646 => "Rwanda", +-99 => "W. Sahara", +682 => "Saudi Arabia", +729 => "Sudan", +728 => "S. Sudan", +686 => "Senegal", +702 => "Singapore", +239 => "S. Geo. and S. Sandw. Is.", +654 => "Saint Helena", +90 => "Solomon Is.", +694 => "Sierra Leone", +222 => "El Salvador", +674 => "San Marino", +-99 => "Somaliland", +706 => "Somalia", +666 => "St. Pierre and Miquelon", +688 => "Serbia", +678 => "Sao Tome and Principe", +740 => "Suriname", +703 => "Slovakia", +705 => "Slovenia", +752 => "Sweden", +748 => "Swaziland", +534 => "Sint Maarten", +690 => "Seychelles", +760 => "Syria", +796 => "Turks and Caicos Is.", +148 => "Chad", +768 => "Togo", +764 => "Thailand", +762 => "Tajikistan", +795 => "Turkmenistan", +626 => "East Timor", +776 => "Tonga", +780 => "Trinidad and Tobago", +788 => "Tunisia", +792 => "Turkey", +798 => "Tuvalu", +158 => "Taiwan", +834 => "Tanzania", +800 => "Uganda", +804 => "Ukraine", +581 => "U.S. Minor Outlying Is.", +858 => "Uruguay", +840 => "United States", +-99 => "Guantanamo Bay USNB", +860 => "Uzbekistan", +336 => "Vatican", +670 => "St. Vin. and Gren.", +862 => "Venezuela", +92 => "British Virgin Is.", +850 => "U.S. Virgin Is.", +704 => "Vietnam", +548 => "Vanuatu", +275 => "West Bank", +876 => "Wallis and Futuna", +-99 => "Akrotiri", +882 => "Samoa", +887 => "Yemen", +710 => "South Africa", +894 => "Zambia", +716 => "Zimbabwe", +typemax(Int32) => "MISSING" +) diff --git a/src/Cubes/Axes.jl b/src/Cubes/Axes.jl index 22e8c2e3..d8d86f51 100644 --- a/src/Cubes/Axes.jl +++ b/src/Cubes/Axes.jl @@ -101,7 +101,6 @@ CategoricalAxis(s::AbstractString,v::Vector)=CategoricalAxis(Symbol(s),v) @defineCatAxis Variable String @defineCatAxis SpatialPoint Tuple{Float64,Float64} -@defineCatAxis Country String @defineCatAxis TimeScale String @defineCatAxis Quantile Float64 diff --git a/src/Cubes/Cubes.jl b/src/Cubes/Cubes.jl index 542b8a31..9386c748 100644 --- a/src/Cubes/Cubes.jl +++ b/src/Cubes/Cubes.jl @@ -79,9 +79,10 @@ type CubeMem{T,N} <: AbstractCubeMem{T,N} axes::Vector{CubeAxis} data::Array{T,N} mask::Array{UInt8,N} + properties::Dict{String} end - +CubeMem(axes::Vector{CubeAxis},data,mask) = CubeMem(axes,data,mask,Dict{String,Any}()) Base.permutedims(c::CubeMem,p)=CubeMem(c.axes[collect(p)],permutedims(c.data,p),permutedims(c.mask,p)) axes(c::CubeMem)=c.axes diff --git a/src/Cubes/TempCubes.jl b/src/Cubes/TempCubes.jl index 9bfa0d85..00abb75b 100644 --- a/src/Cubes/TempCubes.jl +++ b/src/Cubes/TempCubes.jl @@ -29,6 +29,7 @@ type TempCube{T,N} <: AbstractTempCube{T,N} folder::String block_size::CartesianIndex{N} persist::Bool + properties::Dict{String} end "This defines a perumtation of a temporary datacube, as a result from perumtedims on a TempCube" type TempCubePerm{T,N} <: AbstractTempCube{T,N} @@ -36,6 +37,7 @@ type TempCubePerm{T,N} <: AbstractTempCube{T,N} folder::String block_size::CartesianIndex{N} perm::NTuple{N,Int} + properties::Dict{String} end @@ -69,7 +71,7 @@ containing the dimensions of each sub-file the cube is split into. - `persist=true` shall the cubes disk representation be kept or deleted when the cube gets out of scope - `overwrite=false` shall the data in an existing folder be overwritten """ -function TempCube{N}(axlist,block_size::CartesianIndex{N};folder=mktempdir(),T=Float32,persist::Bool=true,overwrite::Bool=false) +function TempCube{N}(axlist,block_size::CartesianIndex{N};folder=mktempdir(),T=Float32,persist::Bool=true,overwrite::Bool=false,properties=Dict{String,Any}()) isdir(folder) || mkpath(folder) if !isempty(readdir(folder)) if overwrite @@ -86,29 +88,30 @@ function TempCube{N}(axlist,block_size::CartesianIndex{N};folder=mktempdir(),T=F ssmall=map(div,s,block_size.I) for ii in CartesianRange(totuple(ssmall)) istart = (CItimes((ii-CartesianIndex{N}()),block_size))+CartesianIndex{N}() - ncdims = NcDim[NcDim(axlist[i],istart[i],block_size[i]) for i=1:N] + ncdims = N>0 ? NcDim[NcDim(axlist[i],istart[i],block_size[i]) for i=1:N] : [NcDim("DummyDim",1)] vars = NcVar[NcVar("cube",ncdims,t=T),NcVar("mask",ncdims,t=UInt8)] nc = NetCDF.create(joinpath(folder,tofilename(ii)),vars) # NetCDF.putvar(nc["cube"],fill(iniVal,block_size)) # NetCDF.putvar(nc["mask"],fill(MISSING,block_size)) NetCDF.close(nc) end - save(joinpath(folder,"axinfo.jld"),"axlist",axlist) - ntc=TempCube{T,N}(axlist,folder,block_size,persist) + save(joinpath(folder,"axinfo.jld"),"axlist",axlist,"properties",properties) + ntc=TempCube{T,N}(axlist,folder,block_size,persist,properties) finalizer(ntc,cleanTempCube) return ntc end -TempCube(axlist,block_size::Tuple;kwargs...)=TempCube(axlist,CartesianIndex(block_size)) +TempCube(axlist,block_size::Tuple;kwargs...)=TempCube(axlist,CartesianIndex(block_size);kwargs...) function openTempCube(folder;persist=true) axlist=load(joinpath(folder,"axinfo.jld"),"axlist") + properties=load(joinpath(folder,"axinfo.jld"),"properties") N=length(axlist) v=NetCDF.open(joinpath(folder,tofilename(CartesianIndex{N}())),"cube") T=eltype(v) - block_size=CartesianIndex(size(v)) + block_size= N==0 ? CartesianIndex(()) : CartesianIndex(size(v)) ncclose(joinpath(folder,tofilename(CartesianIndex{N}()))) - return TempCube{T,N}(axlist,folder,block_size,persist) + return TempCube{T,N}(axlist,folder,block_size,persist,properties) end function readCubeData{T}(y::AbstractTempCube{T}) @@ -149,7 +152,7 @@ function _read{N}(y::TempCube,thedata::NTuple{2},r::CartesianRange{CartesianInde return nothing end -Base.permutedims{T,N}(c::TempCube{T,N},perm)=TempCubePerm{T,N}(c.axes,c.folder,c.block_size,perm) +Base.permutedims{T,N}(c::TempCube{T,N},perm)=TempCubePerm{T,N}(c.axes,c.folder,c.block_size,perm,c.properties) #Method for reading cubes that get transposed #Per means fileOrder -> MemoryOrder function _read{N}(y::TempCubePerm,thedata::NTuple{2},r::CartesianRange{CartesianIndex{N}}) diff --git a/src/DAT/DAT.jl b/src/DAT/DAT.jl index 9060ea07..43825ac4 100644 --- a/src/DAT/DAT.jl +++ b/src/DAT/DAT.jl @@ -357,7 +357,7 @@ function generateOutCube(dc::DATConfig,gfun,T1,i) @inbounds for ii in eachindex(outar) outar[ii]=gfun(T1) end - dc.outcubes[i] = Cubes.CubeMem{T,length(newsize)}(dc.axlistOut[i], outar,zeros(UInt8,newsize...)) + dc.outcubes[i] = Cubes.CubeMem(dc.axlistOut[i], outar,zeros(UInt8,newsize...)) else error("return cube type not defined") end diff --git a/src/Plot/Plot.jl b/src/Plot/Plot.jl index 00e4c7c8..cfe20949 100644 --- a/src/Plot/Plot.jl +++ b/src/Plot/Plot.jl @@ -16,6 +16,7 @@ import Base.Cartesian: @ntuple,@nexprs import Patchwork.load_js_runtime import Measures import Compose +import DataStructures: OrderedDict #import Patchwork.load_js_runtime ga=[] @@ -197,10 +198,20 @@ function val2col(x,m,colorm,mi,ma,misscol,oceancol) end end +function val2col(x,m,colorm::Dict,misscol,oceancol) + if !isnan(x) && m==VALID || m==FILLED + return get(colorm,x,misscol) + elseif (m & OCEAN)==OCEAN + return oceancol + else + return misscol + end +end + import PlotUtils.optimize_ticks import PlotUtils.cgrad import Colors.@colorant_str -import Compose: rectangle, text, line, compose, context, stroke, svgattribute, bitmap, HCenter, VBottom +import Compose: rectangle, text, line, compose, context, stroke, svgattribute, bitmap, HCenter, VBottom, HRight, VCenter const namedcolms=Dict( :viridis=>[cgrad(:viridis)[ix] for ix in linspace(0,1,100)], :magma=>[cgrad(:magma)[ix] for ix in linspace(0,1,100)], @@ -227,8 +238,7 @@ Map plotting tool for cube objects, can be called on any type of cube data If a dimension is neither longitude or latitude and is not fixed through an additional keyword, a slider or dropdown menu will appear to select the axis value. """ function plotMAP{T}(cube::CubeAPI.AbstractCubeData{T};dmin=zero(T),dmax=zero(T), - colorm=:inferno,oceancol=colorant"darkblue",misscol=colorant"gray",symmetric=false, - labels=nothing,kwargs...) + colorm=:inferno,oceancol=colorant"darkblue",misscol=colorant"gray",symmetric=false,kwargs...) isa(colorm,Symbol) && (colorm=namedcolms[colorm]) dmin,dmax=typed_dminmax(T,dmin,dmax) axlist=axes(cube) @@ -249,14 +259,22 @@ function plotMAP{T}(cube::CubeAPI.AbstractCubeData{T};dmin=zero(T),dmax=zero(T), subcubedims[2]=nlat sliceargs=Any[:(1:$nlon),:(1:$nlat)] fixedvarsEx=quote end - if labels!=nothing + if haskey(cube.properties,"labels") + labels = cube.properties["labels"] push!(fixedvarsEx.args,:(labels=$labels)) - colorm=distinguishable_colors(length(labels)+2,[misscol,oceancol])[3:end] - legex = :(getlegend(colorm,legheight,labels)) - mimaex = :(mi=1;ma=$(length(labels))) + _colorm=distinguishable_colors(length(labels)+2,[misscol,oceancol])[3:end] + colorm=Dict(k=>_colorm[i] for (i,k) in enumerate(keys(labels))) + colorm2=Dict(k=>_colorm[i] for (i,k) in enumerate(values(labels))) + rgbarEx = :(rgbar=getRGBAR(a,m,colorm,misscol,oceancol,nx,ny)) + legEx = :(getlegend(colorm2,legwidth)) + mimaex = :(colorm2=$colorm2) + legPosEx = :(legPos=:right) else - legex = :(getlegend(mi,ma,colorm,legheight)) - mimaex = dmin==dmax ? :((mi,ma)=getMinMax(a,m,symmetric=$symmetric)) : :(mi=$(dmin);ma=$(dmax)) + labels=nothing + mimaex = dmin==dmax ? :((mi,ma)=getMinMax(a,m,symmetric=$symmetric)) : :(mi=$(dmin);ma=$(dmax)) + rgbarEx = :(rgbar=getRGBAR(a,m,colorm,convert($T,mi),convert($T,ma),misscol,oceancol,nx,ny)) + legEx = :(getlegend(mi,ma,colorm,legheight)) + legPosEx = :(legPos=:bottom) end fixedAxes=CubeAxis[] for (sy,val) in kwargs @@ -299,12 +317,14 @@ function plotMAP{T}(cube::CubeAPI.AbstractCubeData{T};dmin=zero(T),dmax=zero(T), colorm=$colorm oceancol=$oceancol misscol=$misscol - rgbar=getRGBAR(a,m,colorm,convert($T,mi),convert($T,ma),misscol,oceancol,nx,ny) + $rgbarEx + $legPosEx pngbuf=IOBuffer() show(pngbuf,"image/png",Image(rgbar,Dict("spatialorder"=>["x","y"]))) - legheight=max(0.1*Measures.h,1.6Measures.cm) - themap=obj=compose(context(0,0,1,1Measures.h-legheight),bitmap("image/png",pngbuf.data,0,0,1,1)) - theleg=$legex + legheight=legPos==:bottom ? max(0.1*Measures.h,1.6Measures.cm) : 0Measures.h + legwidth =legPos==:right ? max(0.2*Measures.w,3.2Measures.cm) : 0Measures.w + themap=compose(context(0,0,1Measures.w-legwidth,1Measures.h-legheight),bitmap("image/png",pngbuf.data,0,0,1,1)) + theleg=$legEx compose(context(),themap,theleg) end lambda = Expr(:(->), Expr(:tuple, argvars...),plotfun) @@ -318,6 +338,7 @@ function plotMAP{T}(cube::CubeAPI.AbstractCubeData{T};dmin=zero(T),dmax=zero(T), myfun() end @noinline getRGBAR(a,m,colorm,mi,ma,misscol,oceancol,nx,ny)=RGB{U8}[val2col(a[i,j],m[i,j],colorm,mi,ma,misscol,oceancol) for i=1:nx,j=1:ny] +@noinline getRGBAR(a,m,colorm::Dict,misscol,oceancol,nx,ny)=RGB{U8}[val2col(a[i,j],m[i,j],colorm,misscol,oceancol) for i=1:nx,j=1:ny] import Showoff.showoff function getlegend(xmin,xmax,colm,legheight) @@ -333,6 +354,20 @@ function getlegend(xmin,xmax,colm,legheight) compose(context(0,1Measures.h-legheight,1,legheight),bar,tlabels,dlines) end +function getlegend(colm,width) + texth1=Compose.max_text_extents(Compose.default_font_family, Compose.default_font_size, first(keys(colm)))[2] + texth=texth1*1.05*length(colm) + yoffs=(Measures.h-texth)/2 + yl=texth + ncol=length(colm) + tpos=[(i-0.5)/ncol for i=1:ncol] + r=Compose.circle([0.5],[(i-0.5)/ncol for i in 1:ncol],[max(1/(ncol-1),0.5)]) + f=fill(collect(values(colm))) + bar=compose(context(0.9,yoffs,0.1,yl),r,f,stroke(nothing),svgattribute("shape-rendering","crispEdges")) + tlabels=compose(context(0,yoffs,0.85,yl),Compose.text([1],tpos,collect(keys(colm)),[HRight()],[VCenter()])) + compose(context(1Measures.w-width,0,width,1),bar,tlabels) + end + """ `plotScatter(cube::AbstractCubeData; vsaxis=VariableAxis, alongaxis=0, group=0, xaxis=0, yaxis=0, kwargs...)` @@ -504,16 +539,5 @@ function plotScatter{T}(cube::AbstractCubeData{T};group=0,vsaxis=VariableAxis,xa display(myfun(cube)) end -function getlegend(colm,legheight,labels) - xoffs=0.05 - xl=1-2xoffs - tlabs=labels - tpos=[(i-0.5)/length(tlabs) for i=1:length(tlabs)] - r=rectangle([(i-1)/length(colm) for i in 1:length(colm)],[0],[1/(length(colm))],[1]) - f=fill([colm[div((i-1)*length(colm),length(colm))+1] for i=1:length(colm)]) - bar=compose(context(xoffs,0.35,xl,0.55),r,f,stroke(nothing),svgattribute("shape-rendering","crispEdges")) - tlabels=compose(context(xoffs,0,xl,0.2),Compose.text(tpos,[1],tlabs,[HCenter()],[VBottom()])) - compose(context(0,1Measures.h-legheight,1,legheight),bar,tlabels) -end end diff --git a/src/Proc/OnlineStats.jl b/src/Proc/OnlineStats.jl index 076d7658..d02beac0 100644 --- a/src/Proc/OnlineStats.jl +++ b/src/Proc/OnlineStats.jl @@ -15,8 +15,8 @@ function DATfitOnline{T<:OnlineStat{OnlineStats.ScalarInput}}(xout::AbstractArra end function DATfitOnline{T<:OnlineStat{OnlineStats.ScalarInput}}(xout::AbstractArray{T},maskout,xin,maskin,splitmask,msplitmask,cfun) - for (mi,xi,si) in zip(maskin,xin,splitmask) - (mi & MISSING)==VALID && fit!(xout[cfun(si)],xi) + for (mi,xi,si,m2) in zip(maskin,xin,splitmask,msplitmask) + ((mi | m2) & MISSING)==VALID && fit!(xout[cfun(si)],xi) end end @@ -54,7 +54,7 @@ function DATfitOnline{T<:OnlineStats.VectorInput,U}(xout::AbstractArray{T},masko end function finalizeOnlineCube{T<:OnlineStat,N}(c::CubeMem{T,N}) - CubeMem{Float32,N}(c.axes,map(OnlineStats.value,c.data),c.mask) + CubeMem(c.axes,map(OnlineStats.value,c.data),c.mask) end function mapCube{T<:OnlineStat}(f::Type{T},cdata::AbstractCubeData;by=CubeAxis[],max_cache=1e7,cfun=identity,outAxis=nothing,MDAxis=nothing,kwargs...) @@ -67,7 +67,17 @@ function mapCube{T<:OnlineStat}(f::Type{T},cdata::AbstractCubeData;by=CubeAxis[] end bycubes=filter(i->!in(i,inaxtypes),collect(by)) if length(bycubes)==1 - outAxis==nothing && error("You have to specify an output axis") + if outAxis==nothing + if haskey(bycubes[1].properties,"labels") + idict=bycubes[1].properties["labels"] + axname=get(bycubes[1].properties,"name","Label") + outAxis=CategoricalAxis(axname,collect(values(idict))) + const convertdict=Dict(k=>i for (i,k) in enumerate(keys(idict))) + cfun=x->convertdict[x] + else + error("You have to specify an output axis") + end + end indata=(cdata,bycubes[1]) isa(outAxis,DataType) && (outAxis=outAxis()) outdims=(outAxis,) diff --git a/test/access.jl b/test/access.jl index 72f57608..cc1dc3b6 100644 --- a/test/access.jl +++ b/test/access.jl @@ -49,6 +49,14 @@ llcube = extractLonLats(data1,ll) @test llcube.data[2,:]==data1.data[3,1,:] @test llcube.data[3,:]==data1.data[4,1,:] + +#Test access datacube by region +d3 = getCubeData(c,variable="gross_primary_productivity",region="Cambodia",time=DateTime("2005-01-01")) +@test d3.axes==[LonAxis(102.25:0.25:107.25),LatAxis(14.75:-0.25:10.75)] + + + + #Test saving cubes dire=mktempdir() CABLABdir(dire) diff --git a/test/online.jl b/test/online.jl index ffe252de..81807bf3 100644 --- a/test/online.jl +++ b/test/online.jl @@ -16,26 +16,26 @@ d2=readCubeData(d) @test_approx_eq mean(d2.data,(1,3))[:] oo2.data srand(1) -mask=CubeMem(d2.axes,rand(1:10,size(d2.data)),zeros(UInt8,size(d2.data))) -mask2=CubeMem(d2.axes[1:2],rand(1:10,4,4),zeros(UInt8,4,4)) +mask=CubeMem(d2.axes,rand(1:10,size(d2.data)),zeros(UInt8,size(d2.data)),Dict("labels"=>Dict(i=>string(i) for i=1:10))) +mask2=CubeMem(d2.axes[1:2],rand(1:10,4,4),zeros(UInt8,4,4),Dict("labels"=>Dict(i=>string(i) for i=1:10))) -outAx=CategoricalAxis("LandClass",collect(1:10)) +oogrouped = mapCube(Mean,d2,by=(mask,)) +@test isa(oogrouped.axes[1],CategoricalAxis{String,:Label}) -oogrouped = mapCube(Mean,d2,by=(mask,),outAxis=outAx) - -oogrouped2 = mapCube(Mean,d2,by=(mask2,),outAxis=outAx) +oogrouped2 = mapCube(Mean,d2,by=(mask2,)) for k=1:10 - @test_approx_eq oogrouped.data[k] mean(d2.data[mask.data.==k]) + know = findfirst(j->j==string(k),oogrouped2.axes[1].values) + @test_approx_eq oogrouped.data[know] mean(d2.data[mask.data.==k]) i=find(mask2.data.==k) - if length(i) >0 + if length(i) > 0 i1,i2=ind2sub((4,4),i) dhelp=Float32[] for (j1,j2) in zip(i1,i2) append!(dhelp,d2.data[j1,j2,:]) end - @test_approx_eq mean(dhelp) oogrouped2.data[k] + @test_approx_eq mean(dhelp) oogrouped2.data[know] end end