@@ -7,78 +7,99 @@ using FillArrays
77# maskfolder = "data/forestcompressed"
88
99function meanvote (orbits, significance_thresh= - 1.28 )
10- s,n = 0.0 ,0
10+ s, n = 0.0 , 0
1111 for i in eachindex (orbits)
1212 if orbits[i] != 0 && ! isnan (orbits[i])
1313 s += orbits[i]
1414 n += 1
1515 end
1616 end
17- m = s/ n
18- m < significance_thresh ? 1 : 0
17+ m = s / n
18+ return m < significance_thresh ? 1 : 0
1919end
2020
21- function filtersmallcomps! (xout,xin_unfiltered,forestmask,comborbits,connsize;dims= :,threaded= false )
22- xin = broadcast (xin_unfiltered,forestmask) do x,m
23- ismissing (m) ? zero (x) : x* m
21+ function filtersmallcomps! (
22+ xout, xin_unfiltered, forestmask, comborbits, connsize; dims= :, threaded= false
23+ )
24+ xin = broadcast (xin_unfiltered, forestmask) do x, m
25+ ismissing (m) ? zero (x) : x * m
2426 end
25- x = similar (Array{Float64},(axes (xin,1 ),axes (xin,2 ),Base. OneTo (1 )))
26- for j in axes (x,2 ), i in axes (x,1 )
27- x[i,j, 1 ] = comborbits (view (xin,i,j, :))
27+ x = similar (Array{Float64}, (axes (xin, 1 ), axes (xin, 2 ), Base. OneTo (1 )))
28+ for j in axes (x, 2 ), i in axes (x, 1 )
29+ x[i, j, 1 ] = comborbits (view (xin, i, j, :))
2830 end
2931 lc = label_components (x)
3032 c = counter (lc)
3133 for ix in eachindex (xout)
3234 v = lc[ix]
33- if v== 0 || c[v] < connsize
35+ if v == 0 || c[v] < connsize
3436 xout[ix] = 0
3537 else
3638 xout[ix] = 1
3739 end
3840 end
3941end
4042
41- function postprocess (a,target_array:: YAXArray ,forestmask:: YAXArray , orbitcombine= meanvote;minsize= 30 ,max_cache= 5e8 )
42- nx,ny,nz = size (a)
43- windowsx = DAE. MovingWindow (1 - minsize,1 ,2 * minsize + 1 ,nx,(1 ,nx))
44- windowsy = DAE. MovingWindow (1 - minsize,1 ,2 * minsize + 1 ,ny,(1 ,ny))
43+ function postprocess (
44+ a:: YAXArray ,
45+ target_array:: YAXArray ,
46+ forestmask:: YAXArray ,
47+ orbitcombine= meanvote;
48+ minsize= 30 ,
49+ max_cache= 5e8 ,
50+ )
51+ nx, ny, nz = size (a)
52+ windowsx = DAE. MovingWindow (1 - minsize, 1 , 2 * minsize + 1 , nx, (1 , nx))
53+ windowsy = DAE. MovingWindow (1 - minsize, 1 , 2 * minsize + 1 , ny, (1 , ny))
4554 windowsz = [1 : nz]
46- inar1 = DAE. InputArray (a. data,windows= (windowsx,windowsy,windowsz));
47- inar2 = DAE. InputArray (forestmask. data, windows= (windowsx,windowsy));
48- inars = (inar1,inar2)
49- outchunks = (target_array. chunks. chunks... ,DAE. RegularChunks (1 ,0 ,1 ))
50- outars = DAE. create_outwindows ((nx,ny,1 ),chunks= outchunks);
51- uf = DAE. create_userfunction (filtersmallcomps!,UInt8,is_blockfunction= true ,is_mutating= true ,args= (orbitcombine,minsize))
52- op = DAE. GMDWop (inars,(outars,),uf)
53- plan = DAE. optimize_loopranges (op,max_cache)
54- runner= DAE. LocalRunner (op,plan,(reshape (target_array. data,(nx,ny,1 )),))
55+ inar1 = DAE. InputArray (a. data; windows= (windowsx, windowsy, windowsz))
56+ inar2 = DAE. InputArray (forestmask. data; windows= (windowsx, windowsy))
57+ inars = (inar1, inar2)
58+ outchunks = (target_array. chunks. chunks... , DAE. RegularChunks (1 , 0 , 1 ))
59+ outars = DAE. create_outwindows ((nx, ny, 1 ); chunks= outchunks)
60+ uf = DAE. create_userfunction (
61+ filtersmallcomps!,
62+ UInt8;
63+ is_blockfunction= true ,
64+ is_mutating= true ,
65+ args= (orbitcombine, minsize),
66+ )
67+ op = DAE. GMDWop (inars, (outars,), uf)
68+ plan = DAE. optimize_loopranges (op, max_cache)
69+ runner = DAE. LocalRunner (op, plan, (reshape (target_array. data, (nx, ny, 1 )),))
5570 run (runner)
56- target_array
71+ return target_array
5772end
5873
5974function postprocess (tile:: AbstractString , indir, outpath, maskfolder)
6075 if isdir (outpath)
61- return
76+ return nothing
6277 end
6378
6479 forpath = only (glob (" *$tile *" , maskfolder))
6580 @show forpath
66- allfiles = readdir (indir, join= true )
81+ allfiles = readdir (indir; join= true )
6782 @show allfiles
68- orbitfiles = filter (x-> occursin (tile,string (x)), allfiles)
83+ orbitfiles = filter (x -> occursin (tile, string (x)), allfiles)
6984 @show orbitfiles
70- orbits = filter (x-> occursin (" .zarr" , string (x)), orbitfiles)
85+ orbits = filter (x -> occursin (" .zarr" , string (x)), orbitfiles)
7186 @show orbits
72- orbitname = map (o-> split (basename (o),' _' )[4 ],orbits)
87+ orbitname = map (o -> split (basename (o), ' _' )[4 ], orbits)
7388 d = DD. format (Dim {:orbit} (orbitname))
74- files = DD. DimArray (orbits,d)
89+ files = DD. DimArray (orbits, d)
7590 ds = open_mfdataset (string .(files))
7691 nx, ny = size (ds. layer)
77- outds_skeleton = Dataset (;defo= YAXArray ((ds. X,ds. Y),Fill (UInt8 (0 ),nx,ny),chunks= DAE. GridChunks ((nx,ny),(256 ,256 ))))
78- dsout = savedataset (outds_skeleton,path= outpath,skeleton= true ,overwrite= true )
92+ outds_skeleton = Dataset (;
93+ defo= YAXArray (
94+ (ds. X, ds. Y),
95+ Fill (UInt8 (0 ), nx, ny);
96+ chunks= DAE. GridChunks ((nx, ny), (256 , 256 )),
97+ ),
98+ )
99+ dsout = savedataset (outds_skeleton; path= outpath, skeleton= true , overwrite= true )
79100 forest = Cube (forpath)
80101 # masked = map(*, ds.layer, setchunks(forest,DiskArrays.eachchunk(ds.layer.chunks)))
81- @time postprocess (ds. layer,dsout. defo, forest)
102+ @time postprocess (ds. layer, dsout. defo, forest)
82103 @time PyramidScheme. buildpyramids (outpath)
83104end
84105
0 commit comments