CosmOz Data Processing
Niranjan Wimalathunge
2023-08-16
Introduction
A cosmic-ray soil moisture probe is designed to passively and non-invasively track water content within the uppermost 50cm of soil or the equivalent depth of water – even extending to 20cm of snow. Employing the cosmic-ray methodology, this technology encompasses an averaging range extending laterally up to roughly 300m at sea level, offering an unprecedented scale of observation. The mobile version of it, the rover, acquires neutron count data across the surface at one-minute intervals. However, it is essential to process these neutron counts to eliminate environmental influences and to calibrate them against field soil moisture measurements. The calibration parameter N0 (the neutron intensity above the dry soil) plays a pivotal role in this calibration process.
Libraries
library(raster)
library(readxl)
library(epiR)
library(sp)
library(rgdal)
library(dplyr)
library(sf)
library(automap)
library(mapview)Area of interest (BARS)
## OGR data source with driver: ESRI Shapefile
## Source: "X:\PRJ-SoilWaterNow\data\Aus\Farms&sites\Bars\BARS_boundary.shp", layer: "BARS_boundary"
## with 1 features
## It has 1 fieldsGet data
Use the link to download the datasets and change the file paths accordingly Survey date: 2022-03-28
Fixed variables when process cosmos neutron count
Pref<-1000.71 # (mb)Reference Pressure:Atmospheric Standards formula (en.wikipedia.org/wiki/Standard_atmospheric_pressure)
x1<-1020.5 # (g cm-2) Atmospheric depth:Reference pressure divided by gravity
beta<-0.0076 # (mb-1)Atmospheric attentuation coeffcient:Desilets and Zreda (2003) equation 11
SLref<-1013.25# (mb)Sea level reference pressure:Atmospheric Standards formula (en.wikipedia.org/wiki/Standard_atmospheric_pressure)
x2<-1033.2 # (g cm-2) Reference pressure divided by gravity
beta2<-0.0077 # (mb-1)Desilets and Zreda (2003) equation 11
ScaleEle<-0.91 # Desilets and Zreda (2006) equation 7 modifed and provided by Darin Desilets
ScaleLat<-0.93 # Desilets and Zreda (2006) equation 6
SclaleLat1<-1.069669 # Latitude Scale (as multiplier)Cutoff Rigidity
Cutoff Rigidity Calculator “http://cosmos.hwr.arizona.edu/Util/rigidity.php” For a given lat/lon, this utility calculates the cutoff rigidity (GV) of incoming primary cosmic rays
Rig=4.24 #Cutoff Rigidity at site= 4.24 GVNeutron count intensity from a selected monitoring site
https://www.nmdb.eu/nest/ follow the instruction in the given document
LMKS<-read.csv("X:/PRJ-SoilWaterNow/data/Aus/Farms&sites/Bars/LMKS.csv")[-1]
LMKS$Date<- strptime(LMKS$Date, format="%Y-%m-%d %H:%M:%S")Rover logger field names
## [1] "RecordNum" "Date" "PTB110_mb" "P1_mb" "T1_C"
## [6] "RH1" "T_CS215" "RH_CS215" "Vbat" "N1Cts"
## [11] "N2Cts" "N3Cts" "N4Cts" "N5Cts" "N6Cts"
## [16] "N7Cts" "N8Cts" "N9Cts" "N10Cts" "N11Cts"
## [21] "N12Cts" "N13Cts" "N14Cts" "N15Cts" "N16Cts"
## [26] "N1ETsec" "N2ETsec" "N3ETsec" "N4ETsec" "N5ETsec"
## [31] "N6ETsec" "N7ETsec" "N8ETsec" "N9ETsec" "N10ETsec"
## [36] "N11ETsec" "N12ETsec" "N13ETsec" "N14ETsec" "N15ETsec"
## [41] "N16ETsec" "GpsUTC" "LatDec" "LongDec" "Alt"
## [46] "Qual" "NumSats" "HDOP" "COG" "Speed_kmh"
## [51] "SpeedQuality" "strDate"Rover logger files
Measurements are saved in the logger with the ROV extension. Save them as CSV
cos1<-read.csv("X:/PRJ-SoilWaterNow/data/Aus/Cosmos/RoverData2022/2203282031-2022Bars.csv",header = T);names(cos1)<-names
cos2<-read.csv("X:/PRJ-SoilWaterNow/data/Aus/Cosmos/RoverData2022/2203280451-2022Bars.csv",header = T);names(cos2)<-names
cos<-rbind(cos1,cos2)Combine the Rover logger measurments with the neutron count intensity
cos <- cos[ order(cos$Date , decreasing = T ),]
cos$Date<- strptime(cos$Date, format="%Y/%m/%d %H:%M:%S")
countAverage<-mean(LMKS$countInten)
LMKS$IntenCorr<-LMKS$countInten/countAverage
cos<-inner_join(cos, LMKS, by = "Date")Correct neutron count
The cosmic-ray neutrons is influenced not only by hydrogen in soil water but also by other hydrogen pools. These factors impact the accuracy of neutron counts and needs to be removed.
# fix coordinates
cos<-cos[!(duplicated(cos$LatDec) & cos$LongDec), ]
cos$LatDec<-paste("-",cos$LatDec,sep="")
cos$LatDec<-as.numeric(cos$LatDec)
cos<-as.data.frame(cos)
countMin<-(cos$N1Cts*60/cos$N1ETsec)+(cos$N2Cts*60/cos$N2ETsec)+(cos$N3Cts*60/cos$N3ETsec)+(cos$N4Cts*60/cos$N4ETsec)+(cos$N5Cts*60/cos$N5ETsec)+(cos$N6Cts*60/cos$N6ETsec)+
(cos$N7Cts*60/cos$N7ETsec)+(cos$N8Cts*60/cos$N8ETsec)+(cos$N9Cts*60/cos$N9ETsec)+(cos$N10Cts*60/cos$N10ETsec)+(cos$N11Cts*60/cos$N11ETsec)+(cos$N12Cts*60/cos$N12ETsec)+
(cos$N13Cts*60/cos$N13ETsec)+(cos$N14Cts*60/cos$N14ETsec)+(cos$N15Cts*60/cos$N15ETsec)+(cos$N16Cts*60/cos$N16ETsec)
vp<- (2165*((0.6108*exp((17.27*cos$T_CS215)/(cos$T_CS215+237.3)))*(cos$RH_CS215/100)))/(cos$T_CS215+273.16) #vp(g/m3)
VapCorr<-(1+0.0054*(vp-0))
PressCorr<-exp(0.0077*((cos$PTB110_mb-1013.25)))
RefPress<-1013.25*exp(-1*(9.80665*0.0289644*cos$Alt)/(8.31447*288.15))
AtmosDepth<-RefPress/(9.80655/10)
AtmosAttenCoeff<-0.0054196+0.00022082*Rig+-0.00000051952*Rig^2+(0.0000072062+-0.0000019702*Rig)*AtmosDepth+(-0.0000000098334+0.0000000034201*Rig)*AtmosDepth^2+(0.0000000000049898+-0.0000000000017192*Rig)*AtmosDepth^3
ElevationScale<-exp((AtmosDepth-x2)/(1/AtmosAttenCoeff))
CorrCounts<-((countMin*PressCorr*VapCorr)/cos$IntenCorr)*ElevationScale*SclaleLat1
cos$Counts<-CorrCounts
CorrCounts## [1] 242.6444 243.0246 248.8667 259.0363 234.9191 224.9091 229.2524 257.8756
## [9] 244.2956 267.7671 217.9806 239.2874 231.2779 227.5209 261.7758 262.6125
## [17] 237.8569 253.4428 231.4583 244.4779 233.0098 215.1031 218.8285 190.6336
## [25] 210.3185 209.8917 220.2660 219.3586 224.8905 214.0499 222.1442 221.9809
## [33] 212.0729 246.2297 244.8054 245.0335 247.2659 226.6244 211.2326 247.7784
## [41] 239.8741 241.3586 231.7651 223.2514 238.2838 236.9331 240.7950 226.4803
## [49] 236.6621 254.4578 234.5265 219.0184 229.9210 228.8299 232.5320 226.5915
## [57] 220.5997 230.9271 256.5362 258.3279 204.4936 229.0682 223.3706 212.1348
## [65] 204.1150 224.5611 221.2790 204.9812 204.9311 239.6303 208.5253 242.4989
## [73] 237.8697 201.9617 226.7901 220.3244 224.7634 208.1544 244.0736 223.5549
## [81] 231.0471 219.5841 234.4060 236.3153 231.2885 237.7382 234.0536 234.4147
## [89] 245.6735 224.8537 253.7293 227.4385 238.1026 226.7628 233.8885 223.2720
## [97] 241.8369 241.4231 249.2087 249.7030 232.5111 222.9221 233.1532 236.1486
## [105] 245.3016 217.8351 244.5141 244.0404 225.1068 229.7359 211.4421 225.1741
## [113] 248.3083 229.1944 213.7228 230.8298 228.9159 257.1262 241.9411 244.4132
## [121] 220.7308 223.3857 233.3699 226.9885 233.8415 228.8128 248.8406 236.5597
## [129] 261.8853 234.2482 237.2738 224.1793 238.6209 280.3161 247.8240 264.3734
## [137] 262.5123 262.0806 272.3413 278.8977 241.6779 261.1269 233.2555 241.8094
## [145] 221.9482 226.4325 221.8615 238.5460 227.0401 239.1197 249.5921 270.0589
## [153] 253.7342 269.4585 272.8012 267.0574 259.7871 259.6779 272.7299 268.7378
## [161] 270.0930 264.6934 269.0331 258.5349 256.1372 283.0941 231.8082 248.6447
## [169] 249.5572 244.4285 240.2124 222.2756 263.1067 225.3610 218.9936 241.2473
## [177] 240.7389 224.7806 231.7228 228.6567 219.4278 229.5230 239.6346 242.6566
## [185] 231.0915Map the correct neutron count
e<-extent(Bars)
dem<-readRDS("X:/PRJ-SoilWaterNow/data/Aus/Farms&sites/Bars/dem.rds")
slopeAspect <- raster::terrain(dem, opt = c("slope", "aspect"), unit = "degrees")
S <- SpatialPoints(expand.grid(seq(e@xmin, e@xmax, length=100),
seq(e@ymin, e@ymax, length=100)))
gridded(S) <- TRUE
# dataset to be predicted
newdata<-as.data.frame(raster::extract(stack(dem,slopeAspect),S))
newdata$x<-S@coords[1:10000]
newdata$y<-S@coords[10001:20000]
coordinates(newdata) <- ~x+y
gridded(newdata) <- TRUE
# map neutron count
Ncount<-data.frame(x=cos$LongDec,y=cos$LatDec,count=cos$Counts)
Ncount<-cbind(Ncount,raster::extract(stack(dem,slopeAspect), cbind(Ncount$x,Ncount$y)))
Ncount<-na.omit(Ncount)
coordinates(Ncount) = ~x+y
# universal Kriging
kr.grid <- autoKrige(count ~ dem + slope , Ncount, newdata)## [using universal kriging]r <-raster(kr.grid$krige_output)
raster::crs(r)<-"+proj=longlat +datum=WGS84"
op <- par(mar=rep(2, 4));lim <- par()
Ncount<-raster::crop(r,Bars)
Ncount<-raster::mask(Ncount,Bars)
plot(Ncount,main=expression('Neutron count'));plot(Bars,add=T)
dev.off()## null device
## 1BARS farm data
covariates<-read.csv("X:/PRJ-SoilWaterNow/data/Aus/Farms&sites/Bars/probesSoildata/bars_soil_data_smprobes.csv")
covariates<-covariates[,c(1:3,6:7,11:13)];names(covariates)[1]<-"Probe_SenA"
sites<-read.csv("X:/PRJ-SoilWaterNow/data/Aus/Farms&sites/Bars/probesSoildata/bars_smprobes_id_data.csv")
covariates<-inner_join(sites[,7:9],covariates,by="Probe_SenA");names(covariates)[3:5]<-c("site_ID","Upper_depth","Lower_depth")
# covariates 0-30 cm
covariates1<-filter(covariates,Lower_depth %in% c('10','20','30'))
covariates1<-as.data.frame(covariates1 %>% group_by(site_ID)%>% summarise(across(everything(), mean)))[-c(4,5)]
# add dem, slope to the covariates dataset####
cov<-cbind(covariates1,as.data.frame(raster::extract(stack(dem,slopeAspect),cbind(covariates1$Longitude,covariates1$Latitude))))
coordinates(cov) = ~ Longitude+LatitudeClay
kr.grid <- autoKrige(clay ~ dem + slope , cov, newdata)## [using universal kriging]r <-raster(kr.grid$krige_output)
proj4string(r)<-"+proj=longlat +datum=WGS84"
op <- par(mar=rep(2, 4));lim <- par()
r.crop<-raster::crop(r,Bars)
Clay<-raster::mask(r.crop,Bars)
plot(Clay,main=expression('Clay (0-30cm) %'));plot(Bars,add=T)
dev.off()## null device
## 1Lattice water
# derived from clay use of a PTF
waterLat<-0.00075*Clay+0.0121
plot(waterLat,main=expression('Lattice Water (0-30cm) %'));plot(Bars,add=T)
dev.off()## null device
## 1Carbon
kr.grid <- autoKrige(SOC ~ dem + slope , cov, newdata)## [using universal kriging]r <-raster(kr.grid$krige_output)
proj4string(r)<-"+proj=longlat +datum=WGS84"
op <- par(mar=rep(2, 4));lim <- par()
Carbon<-raster::crop(r,Bars)
Carbon<-raster::mask(Carbon,Bars)
plot(Carbon,main=expression('Carbon (0-30cm) %'));plot(Bars,add=T)
dev.off()## null device
## 1Bulk Density
kr.grid <- autoKrige(bulk_density ~ dem + slope , cov, newdata)## [using universal kriging]r <-raster(kr.grid$krige_output)
proj4string(r)<-"+proj=longlat +datum=WGS84"
op <- par(mar=rep(2, 4));lim <- par()
r.crop<-raster::crop(r,Bars)
BD<-raster::mask(r.crop,Bars)
plot(BD,main=expression('Bulk density (0-30cm) %'));plot(Bars,add=T)
dev.off()## null device
## 1Field soil moisture (0-10cm)
# Field measurements for the calibration
SM<-na.omit(read.csv("X:/PRJ-PAHGISL/Soil/Point/Farms/Boorowa/2022_borrowa_moisture.csv")[,c(1:3,8)])
coords<-read.csv("X:/PRJ-PAHGISL/Soil/Point/Farms/Boorowa/samplecoordinatesBars2023.csv")
SM<-inner_join(SM, coords, by = "Site_ID");names(SM)[4]<-"sm"
SM<-cbind(SM,as.data.frame(raster::extract(stack(dem,slopeAspect),cbind(SM$long,SM$lat))))
coordinates(SM) = ~ long+lat
kr.grid <- autoKrige(sm ~ dem + slope , SM, newdata)## [using universal kriging]r <-raster(kr.grid$krige_output)
proj4string(r)<-"+proj=longlat +datum=WGS84"
op <- par(mar=rep(2, 4));lim <- par()
r.crop<-raster::crop(r,Bars)
SM<-raster::mask(r.crop,Bars)
plot(SM/100,main=expression('Field soil moisture (0-10cm) - Volumetric moisture (cm'^"3"*' cm'^"-3"*")"));plot(Bars,add=T)
dev.off()## null device
## 1Calibration
The corrected neutron counts from the cosmic-ray soil moisture probe are converted to volumetric soil water content using a calibration function
The neutron intensity above the dry soil(N0) is determined by calculating the soil water using the calibration equation farm while changing the N0. Subsequently, the corresponding RMSE (Root Mean Square Error) against the surveyed soil water is calculated. This method is applied to the collected neutron count measurements. The clay, carbon, bulk density, lattice water, and farm soil moisture values are extracted from corresponding interpolated maps.The best N0 is selected, where it gives the lowest RMSE.
# calculate the nnote (N0)
cosmos<-data.frame(x=cos$LongDec,y=cos$LatDec,count=cos$Counts)
cosmos$Carbon<-raster::extract(Carbon,cbind(cosmos$x,cosmos$y))*0.556/100
cosmos$BD<-raster::extract(BD,cbind(cosmos$x,cosmos$y))
cosmos$Clay<-raster::extract(Clay,cbind(cosmos$x,cosmos$y))
cosmos$waterLat<-0.00075*cosmos$Clay+0.0121
cosmos$SM<-(raster::extract(SM,cbind(cosmos$x,cosmos$y))*cosmos$BD)/100
cosmos$CorrCounts<-raster::extract(Ncount,cbind(cosmos$x,cosmos$y))
cosmos<-na.omit(cosmos)
df<-NULL
nnotedf<-NULL
l<-seq(1,1000,1)#first run
for (b in 1:2){
for (a in 1:length(l)){
neutron<-(cosmos$CorrCounts/l[a])-0.372 # calibration fuction
theta<-((0.0808/neutron)-0.115-cosmos$waterLat-(cosmos$Carbon))*cosmos$BD #calibration function
#validation
lins_con <-epi.ccc(theta,cosmos$SM,ci = "z-transform",conf.level = 0.95)
specify_decimal <- function(x, k) format(round(x, k), nsmall=k)
lab <- paste("CCC: ", round(lins_con$rho.c[,1], digits = 2))
rmse<-sqrt(mean((theta-cosmos$SM)^2))
df<-data.frame(l[a],lab,rmse)
nnotedf<-rbind(nnotedf,df)
a=a+1
}
nnote<-nnotedf[nnotedf$rmse==min(nnotedf$rmse),]$l.a
l<-seq(nnote-1,nnote+1,0.001)#second run
}
nnote## [1] 311.522Evaluation of CosmOz measurements
The corresponding lowest RMSE, highest LCCC, and correlation occur at the optimal N0.
neutron<-(cosmos$CorrCounts/nnote)-0.372
cosmos$theta<-((0.0808/neutron)-0.115-cosmos$waterLat-(cosmos$Carbon))*cosmos$BD
lins_con <-epi.ccc(cosmos$theta,cosmos$SM,ci = "z-transform",conf.level = 0.95)
specify_decimal <- function(x, k) format(round(x, k), nsmall=k)
lab <- paste("CCC: ", round(lins_con$rho.c[,1], digits = 2))
rmse<-sqrt(mean((cosmos$theta-cosmos$SM)^2))
plot(cosmos$theta,cosmos$SM,xlim = c(.0, 0.5), ylim = c(.0,0.5),xlab =expression(bold('CosmOz - Volumetric moisture (cm'^"3"*' cm'^"-3"*")")),
ylab =expression(bold('Field - Volumetric moisture (cm'^"3"*' cm'^"-3"*")")), col=rgb(0, 0, 0, 0.7), pch = 16,main="0-10 cm",font=2)
abline(0,1,lwd = 2)
rmse<-sqrt(mean((cosmos$theta-cosmos$SM)^2))
legend(x = "topleft", legend = c("Concordance correlation",lab),bty = "n")
text(0.05,0.42, paste("RMSE =", round(rmse,2)))
text(0.04,0.39, paste("Cor =", round(cor(cosmos$theta,cosmos$SM),2)))
dev.off()## null device
## 1Map the cosmos soil moisture
neutron<-(Ncount/nnote)-0.372
cosmos10<-((0.0808/neutron)-0.115-waterLat-(Carbon*0.556/100))*BD
plot(cosmos10,main=expression('CosmOz (0-10cm) - Volumetric moisture (cm'^"3"*' cm'^"-3"*")"));plot(Bars,add=T)