# title         : rw_HRsoildata.R
# purpose       : Reading and writing of the Croatian soil profile data;
# reference     : The Croatian national soil profile db (2199 profiles);
# producer      : Prepared by T. Hengl
# last update   : In Wageningen, NL, Aug 2012.
# inputs        : "HRSGDB_profiles.txt" all in one table
# outputs       : Horizons and sites data frames;
# remarks 1     : the table data have been published in a book -- Martinovic J, Vrankovic A (1997). Baza tala Republike Hrvatske. I-III. Ministarstvo zaštite okoliša i prostornog planiranja, Zagreb. p. 365., but the data is not publicly available;

download.file("http://globalsoilmap.net/data/HRSGDB_profiles.zip", destfile=paste(getwd(),"HRSGDB_profiles.zip",sep="/"))
unzip("HRSGDB_profiles.zip")

library(aqp)
library(plyr)

## Import the soil profile data HRSDGB:
HRSGDB <- read.delim("HRSGDB_profiles.txt", stringsAsFactors = FALSE) # , stringsAsFactors = FALSE
str(HRSGDB)
HRSGDB$corrX <- ifelse(!is.na(HRSGDB$corrX), HRSGDB$corrX, HRSGDB$Cro16.30_X)
HRSGDB$corrY <- ifelse(!is.na(HRSGDB$corrY), HRSGDB$corrY, HRSGDB$Cro16.30_Y)
HRSGDB <- HRSGDB[!is.na(HRSGDB$corrX)&!is.na(HRSGDB$corrY),]
coordinates(HRSGDB) <- ~ corrX + corrY
HRSGDB$Cro16.30_X <- NULL
HRSGDB$Cro16.30_Y <- NULL
proj4string(HRSGDB) <- CRS("+proj=tmerc +lat_0=0 +lon_0=16.5 +k=0.9999 +x_0=2500000 +y_0=0 +ellps=bessel +towgs84=550.499,164.116,475.142,5.80967,2.07902,-11.62386,0.99999445824 +units=m")
HRSGDB.ll <- spTransform(HRSGDB, CRS("+proj=longlat +datum=WGS84"))
site <- data.frame(HRSGDB.ll)

# rename columns:
names(site)[which(names(site)=="corrX")] <- "LONWGS84" # longitude
names(site)[which(names(site)=="corrY")] <- "LATWGS84" # latitude
names(site)[which(names(site)=="FAO_NAME")] <- "TAXNFAO" # old FAO
site$SOURCEID <- as.factor(paste("HR", site$PROF_ID, sep="_"))
site$TIMESTRR <- as.Date(as.character(site$UZORAK), format="%Y")
site$SOURCEDB <- "HRSPDB"

## sort horizons:
# TH: Silt is in 0.02-0.002 mm and not in the 0.05-0.002 system!!
hor1 <- data.frame(SOURCEID=site$SOURCEID, UHDICM=as.numeric(site$HOR1_GOR), LHDICM=as.numeric(site$HOR1_DON), PHIHO5=as.numeric(site$HOR1_PH1), PHIKCL=as.numeric(site$HOR1_PH2), SNDPPT=as.numeric(site$HOR1_MSP)+as.numeric(site$HOR1_MKP), SLTPPT=as.numeric(site$HOR1_MP), CLYPPT=as.numeric(site$HOR1_MG), ORCDRC=as.numeric(site$HOR1_HUM)*10/1.724)
hor2 <- data.frame(SOURCEID=site$SOURCEID, UHDICM=as.numeric(site$HOR2_GOR), LHDICM=as.numeric(site$HOR2_DON), PHIHO5=as.numeric(site$HOR2_PH1), PHIKCL=as.numeric(site$HOR2_PH2), SNDPPT=as.numeric(site$HOR2_MSP)+as.numeric(site$HOR2_MKP), SLTPPT=as.numeric(site$HOR2_MP), CLYPPT=as.numeric(site$HOR2_MG), ORCDRC=as.numeric(site$HOR2_HUM)*10/1.724)
hor3 <- data.frame(SOURCEID=site$SOURCEID, UHDICM=as.numeric(site$HOR3_GOR), LHDICM=as.numeric(site$HOR3_DON), PHIHO5=as.numeric(site$HOR3_PH1), PHIKCL=as.numeric(site$HOR3_PH2), SNDPPT=as.numeric(site$HOR3_MSP)+as.numeric(site$HOR3_MKP), SLTPPT=as.numeric(site$HOR3_MP), CLYPPT=as.numeric(site$HOR3_MG), ORCDRC=as.numeric(site$HOR3_HUM)*10/1.724)
hor4 <- data.frame(SOURCEID=site$SOURCEID, UHDICM=as.numeric(site$HOR4_GOR), LHDICM=as.numeric(site$HOR4_DON), PHIHO5=as.numeric(site$HOR4_PH1), PHIKCL=as.numeric(site$HOR4_PH2), SNDPPT=as.numeric(site$HOR4_MSP)+as.numeric(site$HOR4_MKP), SLTPPT=as.numeric(site$HOR4_MP), CLYPPT=as.numeric(site$HOR4_MG), ORCDRC=as.numeric(site$HOR4_HUM)*10/1.724)
hor5 <- data.frame(SOURCEID=site$SOURCEID, UHDICM=as.numeric(site$HOR5_GOR), LHDICM=as.numeric(site$HOR5_DON), PHIHO5=as.numeric(site$HOR5_PH1), PHIKCL=as.numeric(site$HOR5_PH2), SNDPPT=as.numeric(site$HOR5_MSP)+as.numeric(site$HOR5_MKP), SLTPPT=as.numeric(site$HOR5_MP), CLYPPT=as.numeric(site$HOR5_MG), ORCDRC=as.numeric(site$HOR5_HUM)*10/1.724)
hor6 <- data.frame(SOURCEID=site$SOURCEID, UHDICM=as.numeric(site$HOR6_GOR), LHDICM=as.numeric(site$HOR6_DON), PHIHO5=as.numeric(site$HOR6_PH1), PHIKCL=as.numeric(site$HOR6_PH2), SNDPPT=as.numeric(site$HOR6_MSP)+as.numeric(site$HOR6_MKP), SLTPPT=as.numeric(site$HOR6_MP), CLYPPT=as.numeric(site$HOR6_MG), ORCDRC=as.numeric(site$HOR6_HUM)*10/1.724)
hor7 <- data.frame(SOURCEID=site$SOURCEID, UHDICM=as.numeric(site$HOR7_GOR), LHDICM=as.numeric(site$HOR7_DON), PHIHO5=as.numeric(site$HOR7_PH1), PHIKCL=as.numeric(site$HOR7_PH2), SNDPPT=as.numeric(site$HOR7_MSP)+as.numeric(site$HOR7_MKP), SLTPPT=as.numeric(site$HOR7_MP), CLYPPT=as.numeric(site$HOR7_MG), ORCDRC=as.numeric(site$HOR7_HUM)*10/1.724)
hor8 <- data.frame(SOURCEID=site$SOURCEID, UHDICM=as.numeric(site$HOR8_GOR), LHDICM=as.numeric(site$HOR8_DON), PHIHO5=as.numeric(site$HOR8_PH1), PHIKCL=as.numeric(site$HOR8_PH2), SNDPPT=as.numeric(site$HOR8_MSP)+as.numeric(site$HOR8_MKP), SLTPPT=as.numeric(site$HOR8_MP), CLYPPT=as.numeric(site$HOR8_MG), ORCDRC=as.numeric(site$HOR8_HUM)*10/1.724)

horizon <- do.call(rbind, list(hor1, hor2, hor3, hor4, hor5, hor6, hor7, hor8))

# filter strange values:
summary(horizon$PHIHO5)
horizon$PHIHO5 <- ifelse(horizon$PHIHO5 < 2 | horizon$PHIHO5 > 10, NA, horizon$PHIHO5)
horizon$PHIKCL <- ifelse(horizon$PHIKCL < 2 | horizon$PHIKCL > 10, NA, horizon$PHIKCL)
summary(horizon$ORCDRC)
# subset to complete data:
sel.h <- !is.na(horizon$LHDICM)&!is.na(horizon$UHDICM)
sel.s <- !is.na(site$LONWGS84)&!is.na(site$LATWGS84)

profs.f <- join(site[sel.s,c("SOURCEID", "LONWGS84", "LATWGS84", "TIMESTRR", "TAXNFAO", "SOURCEDB")], horizon[sel.h,], type='inner')
depths(profs.f) <- SOURCEID ~ UHDICM + LHDICM
# extract site data
site(profs.f) <- ~ LONWGS84 + LATWGS84 + TIMESTRR + TAXNFAO + SOURCEDB
## check if there are still some duplicates
profs.f@site$SOURCEID[duplicated(profs.f@site$SOURCEID)]
## spatial duplicates:
sp <- profs.f@site
coordinates(sp) <- ~ LONWGS84 + LATWGS84
sp::zerodist(sp)
## 18 duplicates!

# prepare a SITE and hor table table:
wsp4 <- list(sites=profs.f@site, horizons=profs.f@horizons)
str(wsp4)
save(wsp4, file="../wsp4.rda", compress="xz")

## Visualize all points:
library(plotKML)
shape = "http://maps.google.com/mapfiles/kml/pal2/icon18.png"
kml(HRSGDB.ll, shape = shape, labels = PROF_ID, kmz = TRUE) 

# end of script;