# NOTE: improve structure, reusability and readability!!!
# - first do cat based calculations (how to guess required xylim?)
# - afterwards do all the plotting
<- (x, SensorName, rn = , MyMap = , = ("CE", "wCE", "uCE"),
use.avg = , use.sym = , use.var = , wTDcutoff = , origin = ,
dx = 2, dy = dx, = (0.01, 0.1, 0.5, 0.9), breaks = (x) ((0, (x)), ),
= (-100, 100), = (-100, 100), add = , alpha = 0.3, axs = ("r", "i"),
main = , = 1, fill = , = , bg.col = , fp_only = , addSource = !fp_only,
showMax = , showSensor = !fp_only, dispSname = showSensor, N0 = , lpos = ,
showPerc = , leg.bg.col = (0.9), cpal = , lty = 1, lwd = 1, decPlaces = 0, sigNums = 1,
addWR = , WRpos = 2, WRfrac = 20, WRscale = 1, xy_transform = , transformArgs = ,
useSTRtree = , avoidGEOS = , addSB = , SBpos = 3, StaticMapArgs = ,
showLegend = !fp_only, verbose = ){
if (fill) {
(!("maptools", quietly = )){
('filled polygons are not yet supported without the retired(!) package maptools')
}
(!("rgeos", quietly = )){
('filled polygons are not yet supported without the retired(!) package rgeos')
}
(!("sp", quietly = )){
("package 'sp' is required for filled polygons - please install the package by running install.packages('sp')")
}
}
# Notizen:
# - Sources spaeter hinzufuegen
# - Methode fuer Catalog hinzufuegen
# - arguments besser organisieren
# - main auch bei MyMap
# - xlim/ylim xy_transform WGS84 etc... besser organisieren!!!
# - xy_transform --> coords_2_WGS84 oder so
# - TrueProj als Argument, oder besser: if add=TRUE, par()$usr -> MyMap$size???
sx <- ((x))
((x, "footprint")){
# breaks=rev(brks)
xy_transform <- x$xy_transform
transformArgs <- x$transformArgs
SensorPosition <- x$SensorPosition
WDmean <- x$WDmean
xlimOriginal <- x$xlimOriginal
ylimOriginal <- x$ylimOriginal
xm <- x$xm
ym <- x$ym
if (!(origin)) {
old_ll <- (x = ym[1], y = xm[1])
# transf to origin
new_ll <- (old_ll, Angle = -WDmean)
# shift origin
new_ll <- (new_ll, 2, (origin))
# transf back
new_ll <- (new_ll, Angle = WDmean)
# shift xm & ym
xm <- xm - old_ll[2]] + new_ll[2]]
ym <- ym - old_ll[1]] + new_ll[1]]
# shift xlimOriginal & ylimOriginal
xlimOriginal <- xlimOriginal - origin[1]
ylimOriginal <- ylimOriginal - origin[2]
}
xy <- x$xy
(!add){
("argument 'add' must be TRUE when providing a footprint object")
}
} {
# convert old versions
x <- (x)
(x)
(x)
(((x, "Version"))){
(("Object '", sx[min((sx), 2)], "' has not yet been converted to version 4.2+"))
(x)
}
# checks:
((SensorName)){
("SensorName must be provided")
} {
getThisSensor <- (SensorName[1])
(!(x[Sensor==getThisSensor]))("No results for Sensor: ",getThisSensor," availabe!")
}
([1] %in% ("CE","wCE","uCE"))
((rn)){
getThisRows <- x[Sensor==getThisSensor,(rn)]
} {
getThisRows <- rn
}
getThisSource <- x[rn%in%getThisRows&Sensor==getThisSensor,Source[1]]
index <- x,(rn%in%getThisRows&Sensor==getThisSensor&Source==getThisSource)]
WD <- x[index,WD]
((N0)){
N0 <- x[index[1],N0]
} {
((x[index,N0]<N0))("Choose N0 smaller or equal to original N0!")
}
Ustar <- x[index,Ustar]
Suu <- x[index,sUu]
Svu <- x[index,sVu]
bw <- x[index,bw]
z <- x[index,SensorHeight]
L <- x[index,L]
Zo <- x[index,Zo]
<- (x, 'ModelInput')'Sources']]
if (()) {
addSource <-
} {
Source <- (procSources(),"SourceList")
}
((index)>1){
WDmean <- (360 + (((WD/180*)*Ustar)/(Ustar),((WD/180*)*Ustar)/(Ustar))/*180) %% 360
} {
WDmean <- WD
}
Calc.Sensors <- procSensors((x,"ModelInput")$)$"Calc.Sensors"
SensorPosition <- Calc.Sensors[Calc.Sensors, "Sensor Name"] %chin% SensorName,
("Sensor Name", "x-Coord (m)", "y-Coord (m)", "Sensor ID", "Node")]
SPR <- (SensorPosition,2:3],Angle=-WDmean)
((origin)){
origin <- (SensorPosition,2:3])
}
(!add){
((main))main <- ([1],
"CE"="C-Footprint",
"wCE"="w'C'-Footprint",
"uCE"="u'C'-Footprint")
((MyMap)){
(SensorPosition,2],SensorPosition,3],="n",= + origin[1],= + origin[2],xlab="",ylab="",=,yaxt="n",xaxt="n",main=main,=)
pr <- ()
(!(bg.col))(pr$usr[c(1,2,2,1)],pr$usr[c(3,3,4,4)],=bg.col)
(1,at=atx<-(pr$usr[1:2]))
(2,at=aty<-(pr$usr[3:4]))
} {
(PlotOnStaticMap,((MyMap=MyMap),StaticMapArgs))
add <-
}
}
pr <- ()
(!(MyMap)){
(!("TrueProj" %in% (StaticMapArgs))||!StaticMapArgs$TrueProj){
xy_transformOriginal <- xy_transform
((xy_transform)){
xy_transform <- (x,y,mymap=MyMap,){
LatLon2XY.centered(mymap,y,x)
}
} {
xy_transform <- (x,y,mymap=MyMap,){
xy <- xy_transformOriginal(x,y,)
LatLon2XY.centered(mymap,xy$y,xy$x)
}
}
}
}
(addSource){
(i Source){
pg <- i
(!(xy_transform)){
pg[1:2] <- (xy_transform,((x=pg,1],y=pg,2]),transformArgs))
}
upg <- (pg,3])
(j upg){
ind <- (pg,3]==j)
(pg[ind,1],pg[ind,2],="#006837")
}
}
}
# transform x/y:
(!(xy_transform)){
SensorPosition, 2:3] <- (xy_transform,((x=SensorPosition,2],y=SensorPosition,3]),transformArgs))
}
(showSensor){
sp_dt <- setnames(as.data.table((x,"ModelInput")$,
("Sensor Name", "x-Coord (m)", "y-Coord (m)", "Sensor ID", "Node")]),
("name", "x", "y", "id", "node"))[name == getThisSensor]
# cat("hier stimmt's noch nicht falls google maps!\n")
(dispSname){
(sp_dt, ((id)) > 1]){
sp_dt, {
ind <- (node)
(x[ind], y[ind], lty = 2, cex = 1.5)
(x[ind], y[ind], = 20, cex = 1.5)
((x), (y), ([1], " (", .BY$id, ")"), cex = 0.6, pos = 4)
}, = id]
} {
sp_dt, {
ind <- (node)
(x[ind], y[ind], lty = 2, cex = 1.5)
(x[ind], y[ind], = 20, cex = 1.5)
((x), (y), [1], cex = 0.6, pos = 4)
}, = id]
}
} {
sp_dt, {
ind <- (node)
(x[ind], y[ind], lty = 2, cex = 1.5)
(x[ind], y[ind], = 20, cex = 1.5)
}, = id]
}
}
<- pr$usr[1:2]
<- pr$usr[3:4]
xlimOriginal <-
ylimOriginal <-
(!(xy_transform)){
parscale1 <- (([1],[1]))
parscale2 <- (([2],[2]))
lim1 <- ((-1,-1),(x)((((xy_transform,((x=x[1],y=x[2]),transformArgs))) - ([1],[1]))^2),control=(reltol=($double.eps)/parscale1))$
lim2 <- ((1,1),(x)((((xy_transform,((x=x[1],y=x[2]),transformArgs))) - ([2],[2]))^2),control=(reltol=($double.eps)/parscale2))$
<- (lim1[1],lim2[1])
<- (lim1[2],lim2[2])
}
xylim <- (x=(,()),y=(,each=2))
xylim <- (xylim,Angle=-WDmean)
xb <- ((xylim,1])-10*dx,(xylim,1])+10*dx,dx)
yb <- ((xylim,2])-10*dy,(xylim,2])+10*dy,dy)
xm <- xb-1] - dx/2
ym <- yb-1] - dy/2
xy <- (0,=(xm),=(ym))
# loop over index
(j (index)){
if (verbose) (j,"/",(index),"\n")
<- (x,index[j])
CatNames <- ,Catalog]
Subset_seed <- , seed]
dummy <- ""
Ustardummy <- -9999
# loop over point sensors
(i (CatNames)){
# i <- 1
xytemp <- (0,=(xm),=(ym))
(dummy!=CatNames[i]|x[index[j],Ustar!=Ustardummy]){
dummy <- CatNames[i]
Catalog <- (CatNames[i])
(x[index[j],], Catalog=Catalog)
uvw <- (Catalog)
# subset?
((Catalog,"N0")>N0){
env <- ()
oseed <- env$.Random.seed
(Subset_seed,kind="L'Ecuyer-CMRG")
takeSub <- (((Catalog,"N0"),N0))
if ((oseed)) {
( = ".Random.seed", envir = env)
} {
(".Random.seed", value = oseed, envir = env)
}
attCat <- (Catalog)
indexNew <- 1:N0
(indexNew) <- takeSub
Catalog <- Catalog[Traj_ID %in% takeSub,]
Catalog,":="(Traj_ID=indexNew[as.character(Traj_ID)])]
(ac in ((attCat) %w/o% ("names","row.names",".internal.selfref","uvw0")))setattr(Catalog,ac,attCat[ac]])
uvw <- uvw[takeSub,]
setattr(Catalog,"uvw0",uvw)
HeaderSub <- ((Catalog,"header"),"\n*** Subset of ",N0," Trajectories ***\n\n",sep="")
(HeaderSub) <- ("TDhead","character")
setattr(Catalog,"N0",N0)
setattr(Catalog,"header",HeaderSub)
}
(!(wTDcutoff)) Catalog[wTD<wTDcutoff,wTD:=wTDcutoff]
Catalog,CE:=2/wTD/N0]
(use.avg)Catalog,CE:=(CE)]
([1]!="CE"){
([1]=="wCE"){
w0 <- (Catalog),"w0"]
} {
w0 <- (Catalog),"u0"] - (Ustar, Zo, L, z)
}
Catalog,CE:=w0[Traj_ID]*CE]
}
(use.sym){
Catalog <- Catalog[rep(1:.N,2)]
Catalog[1:(.N/2),y:=-y]
Catalog,CE:=CE/2]
}
}
Ctlg <- (Catalog)
(Ctlg,WD[j]-WDmean)
Ctlg <- Ctlg,":="(x=x+SPR[i,1],y=y+SPR[i,2],CE=CE/(CatNames)/(index))]
x>=xb[1]&x<=xb[length(xb)]&y>=yb[1]&y<=yb[length(yb)]
]
Ctlg,":="(xtag=(x,xb, rightmost.closed = ),
ytag=(y,yb, rightmost.closed = ))]
Cat <- Ctlg,(CE),=.(xtag,ytag)]
xytemp[Cat,(xtag,ytag)]] <- Cat,V1]
(use.var){
xb2l <- Ctlg,(((SPR[i,1]-3*dx/2,(x,SPR[i,1])-10*dx,-2*dx)),(SPR[i,1]+dx/2,(x,SPR[i,1])+10*dx,2*dx))]
yb2l <- Ctlg,(((SPR[i,2]-3*dy/2,(y,SPR[i,2])-10*dy,-2*dy)),(SPR[i,2]+dy/2,(y,SPR[i,2])+10*dy,2*dy))]
xb2r <- xb2l + dx
yb2r <- yb2l + dx
Ctlg,":="(
xtag_ll=(x,xb2l),ytag_ll=(y,yb2l)
,xtag_lr=(x,xb2l),ytag_lr=(y,yb2r)
,xtag_rl=(x,xb2r),ytag_rl=(y,yb2l)
,xtag_rr=(x,xb2r),ytag_rr=(y,yb2r)
)]
Cat_ll <- Ctlg,(CE),=.(xtag_ll,ytag_ll)]
Cat_lr <- Ctlg,(CE),=.(xtag_lr,ytag_lr)]
Cat_rl <- Ctlg,(CE),=.(xtag_rl,ytag_rl)]
Cat_rr <- Ctlg,(CE),=.(xtag_rr,ytag_rr)]
xy_ll <- xy_lr <- xy_rl <- xy_rr <- (0,=(xb2r)-1,=(yb2r)-1)
xy_ll[Cat_ll,(xtag_ll,ytag_ll)]] <- Cat_ll,V1]
xy_lr[Cat_lr,(xtag_lr,ytag_lr)]] <- Cat_lr,V1]
xy_rl[Cat_rl,(xtag_rl,ytag_rl)]] <- Cat_rl,V1]
xy_rr[Cat_rr,(xtag_rr,ytag_rr)]] <- Cat_rr,V1]
xmil <- (xm,xb2l,all.inside = )
ymil <- (ym,yb2l,all.inside = )
xmir <- (xm,xb2r,all.inside = )
ymir <- (ym,yb2r,all.inside = )
xyi_ll <- (expand.grid(x=xmil,y=ymil,KEEP.OUT.ATTRS=))
xyi_rr <- (expand.grid(x=xmir,y=ymir,KEEP.OUT.ATTRS=))
xyi_lr <- (expand.grid(x=xmil,y=ymir,KEEP.OUT.ATTRS=))
xyi_rl <- (expand.grid(x=xmir,y=ymil,KEEP.OUT.ATTRS=))
Mu <- (xy_ll[xyi_ll] + xy_lr[xyi_lr] + xy_rl[xyi_rl] + xy_rr[xyi_rr])/4
Var <- (4*xytemp - Mu)^2
xytemp <- (xytemp*Var + ((Var)-Var)/4*Mu)/((Var))
}
xy <- xy + xytemp/(dx*dy)
}
}
(Catalog,Ctlg,Cat)
xy,(1,(xy))] <- 0
xy[c(1,(xy)),] <- 0
}
brks <- breaks(xy)
cl <- (xm,ym,xy,=brks)
clRot <- (cl,,Angle=WDmean)
clFac <- ((clRot,"[[","level"))
uclFac <- (clFac)
# transform x/y:
(!(xy_transform)){
clRot <- (clRot,(x){
a <- (xy_transform,((x=x$x,y=x$y),transformArgs))
(level=x$level,x=a[1]],y=a[2]])
})
}
if (fill) {
clP <- (clRot,(x)sp::Polygon((x$x,x$y)))
clPs <- (uclFac,(x,y,z,a,b)maptools::checkPolygonsHoles(sp::Polygons(y[z==x],(x)),useSTRtree=a,avoidGEOS=b),y=clP,z=clFac,a=useSTRtree,b=avoidGEOS)
clSp <- sp::SpatialPolygons(clPs,(uclFac))
# Loesung:
mlt <- (axs[1],
"i"=1.08,
1
)
xlim2 <- xlimOriginal - (xlimOriginal)*(1-1/mlt)/2;ylim2 <- ylimOriginal - (ylimOriginal)*(1-1/mlt)/2
c1 <- sp::Polygon(((xlim2[1],xlim2[2],xlim2[2],xlim2[1], xlim2[1]),(ylim2[1],ylim2[1],ylim2[2],ylim2[2],ylim2[1])))
c2 <- sp::Polygons((c1), "sclip")
Pclip <- sp::SpatialPolygons((c2))
SpPclip <- rgeos::gIntersection(clSp, Pclip, byid = )
}
alphachar <- ((alpha*255))
((cpal)){
cpal <- ([1]=="CE") ((brks)) ((brks))
} {
cpal <- cpal((brks))
}
(!add){(h=aty,lty="dotted",="lightgrey");(v=atx,lty="dotted",="lightgrey")}
if (fill) {
fillcol <- (cpal, alphachar)
sp::(SpPclip, = fillcol, border = cpal, lty = lty, lwd = lwd, add = )
()
} {
for (icl clRot) {
ibreak <- (icl$level == brks)
(icl$x, icl$y,
= cpal(ibreak - 1) %% (cpal) + 1],
lty = lty(ibreak - 1) %% (lty) + 1],
lwd = lwd(ibreak - 1) %% (lwd) + 1]
)
}
}
mi <- (xy==(xy),arr.ind=T)
(showMax){
maxli <- ((x=xm[mi[1]],y=ym[mi[2]]),Angle=WD+90,Center=origin)+(SensorPosition,2],SensorPosition,3])
(maxli,=20,cex=0.4,="red")
(maxli[1],maxli[2],"max.",cex=0.6,pos=4,=0.2,="red")
}
(showLegend){
((lpos))lpos <- ("bottomleft","topleft","topright","bottomright")[floor(WDmean/90)+1]
(showPerc){
ltex <- (("%2.",decPlaces,"f%%"),(brks)/(xy)*100)
ltex <- (">",ltex)
if (fill) {
temp <- (lpos,=(" ",(ltex)),="% of max.",text.width=(("00",(decPlaces>0){(".",(("0",decPlaces),collapse=""))}{""}," %")),fill=((cpal),alphachar),border=(cpal),bg=leg.bg.col)
} {
temp <- (lpos, = (" ", (ltex)), = "% of max.",
text.width = (("00",
if (decPlaces > 0) (".", (("0", decPlaces), collapse = "")) "", " %")),
= (cpal), bg = leg.bg.col, lty = (lty), lwd = (lwd))
}
(temp$$left + temp$$w, temp$$y,ltex, pos = 2)
} {
ltex <- ((brks),sigNums)
ltex <- (">",ltex)
if (fill) {
(lpos, = ltex, xjust = 0, fill = ((cpal), alphachar), border = (cpal), bg = leg.bg.col)
} {
(lpos, = ltex, xjust = 0, = (cpal), bg = leg.bg.col, lty = (lty), lwd = (lwd))
}
}
}
(addWR)(WDmean,pos=WRpos,=WRfrac,scF=WRscale)
<- ((MyMap)) 1 MyMap
(addSB)(SBpos,)
wm <- (xy==(xy),arr.ind=)
WDrad <- WDmean/180*
out <- ((
breaks=(brks),
range.z=(xy),
max.index = mi,
xy_transform = xy_transform,
transformArgs = transformArgs,
SensorPosition = SensorPosition,
WDmean = WDmean,
xlimOriginal = xlimOriginal,
ylimOriginal = ylimOriginal,
xm = xm,
ym = ym,
xy = xy), = "footprint"
)
((out))
}
