These functions use the empircally motivated cosmic variance percentage formula provided in Driver & Robotham (2010) Eqn 4.
cosvarsph is a 'best effort' approximation of the comoving box subtended by the specified spherical coordinates using the following conversions:
CoDistLow = cosdistCoDist(z=zmin,H0=70,OmegaM=0.3)
CoDistHigh = cosdistCoDist(z=zmax,H0=70,OmegaM=0.3)
cside=CoDistHigh-CoDistLow
area=skyarea(long = long, lat = lat, inunit = inunit, outunit='deg2')[1]
volume=cosvol(area=area, zmax = zmax, zmin=zmin, H0 = 70, OmegaM = 0.3, inunit='deg2')[1]
aside=cos(mean(lat)*pi/180)*(abs(diff(long))/360)*(CoDistLow+cside/2)
bside=(abs(diff(long))/180)*(CoDistLow+cside/2)
scale=sqrt(volume*1e9/(aside*bside*cside))
aside=aside*scale
bside=bside*scale
return(cosvarcar(aside=aside, bside=bside, cside=cside, subsets=subsets))
cosvararea is a simplifed version of cosvarsph, where the assumption is that aside=bside (so the aspect ratio on the sky is 1:1).