Upload New File

R_scripts/WoSIS_Latest_With_R.R

+## ----setup, include=FALSE---------------------------------------------------------------------------------------------------------
+knitr::opts_chunk$set(echo = TRUE, warnings = FALSE, purl = FALSE)
+
+
+## ----load.package-----------------------------------------------------------------------------------------------------------------
+library(rgdal)          # GDAL access from R
+library(gdalUtils)      # wrappers for GDAL utility programs that could be
+                        #  called from the command line
+library(sp)             # spatial data types, ASDAR structures
+library(sf)             # spatial data types -- Simple Features
+library(ggplot2)        # gpplot graphics
+
+
+
+
+## ---------------------------------------------------------------------------------------------------------------------------------
+drivers <- ogrDrivers()
+# print(drivers)
+ix <- grep("GPKG", drivers$name,  fixed=TRUE)
+drivers[ix,]
+ix <- grep("ESRI", drivers$name,  fixed=TRUE)
+drivers[ix,]
+
+
+## ----specify.wfs.address----------------------------------------------------------------------------------------------------------
+wfs <- "WFS:http://maps.isric.org/mapserv?map=/map/wosis_latest.map"
+
+
+## ----ogrinfo.wfs------------------------------------------------------------------------------------------------------------------
+system.time(
+  layers.info <- ogrinfo(wfs, ro=TRUE, so=TRUE, q=TRUE)
+)
+cat(layers.info, sep = "\n")
+
+
+## ----ogrinfo.site-----------------------------------------------------------------------------------------------------------------
+system.time(
+  profiles.info <-
+    ogrinfo(wfs, layer="ms:wosis_latest_profiles",
+            ro=TRUE, so=TRUE, q=TRUE)
+)
+cat(profiles.info, sep="\n")
+
+
+## ----ogrinfo.eu.profiles.1--------------------------------------------------------------------------------------------------------
+system.time(
+  central.eu.profiles.info <-
+    ogrinfo(wfs, ro=TRUE, so=TRUE, q=FALSE,
+            layer="ms:wosis_latest_profiles",
+            spat=c(6, 48, 8, 50))
+)
+head(central.eu.profiles.info, 8)
+central.eu.profiles.info[41:60]
+
+
+## ----ogrinfo.eu.profiles.2--------------------------------------------------------------------------------------------------------
+ix.f <- grep("Feature Count", central.eu.profiles.info)
+central.eu.profiles.info[ix.f]
+ix.e <- grep("Extent", central.eu.profiles.info)
+central.eu.profiles.info[ix.e]
+ix.g <- grep("GEOGCRS", central.eu.profiles.info)
+cat(paste(central.eu.profiles.info[ix.g:(ix.g+17)], collapse="\n"))
+
+
+## ----show.fields------------------------------------------------------------------------------------------------------------------
+ix.p <- grep("Geometry Column", central.eu.profiles.info)
+n <- length(central.eu.profiles.info)
+central.eu.profiles.info[ix.p+1:n]
+
+
+## ----ogrinfo.india.profiles-------------------------------------------------------------------------------------------------------
+system.time(
+  india.profiles.info <-
+    ogrinfo(wfs, ro=TRUE, so=TRUE, q=FALSE,
+            layer="ms:wosis_latest_profiles",
+            where="country_name='India'")
+)
+
+
+## ----ogrinfo.india.profiles.2-----------------------------------------------------------------------------------------------------
+ix.f <- grep("Feature Count", india.profiles.info)
+india.profiles.info[ix.f]
+ix.e <- grep("Extent", india.profiles.info)
+india.profiles.info[ix.e]
+
+
+## ----ogrinfo.properties, warnings-------------------------------------------------------------------------------------------------
+system.time(
+  property.info <- ogrinfo(wfs, layer="ms:wosis_latest_bdfi33",
+                           ro=TRUE, so=TRUE, q=FALSE)
+)
+cat(property.info, sep="\n")
+ix.f <- grep("Feature Count", property.info)
+property.info[ix.f]
+
+
+## ----ogrinfo.bd.india-------------------------------------------------------------------------------------------------------------
+system.time(
+  bd.india.info <- ogrinfo(wfs, layer="ms:wosis_latest_bdfi33",
+                           where="country_name='India' AND upper_depth > 100",
+                           ro=TRUE, so=TRUE, q=FALSE)
+)
+ix.f <- grep("Feature Count", bd.india.info)
+bd.india.info[ix.f]
+(n.records <- as.numeric(strsplit(bd.india.info[ix.f],
+                                 split=": ", fixed=TRUE)[[1]][2]))
+
+
+## ----setup.local.dir--------------------------------------------------------------------------------------------------------------
+wosis.dir.name <- "./wosis_latest"
+if (!file.exists(wosis.dir.name)) dir.create(wosis.dir.name)
+
+
+## ----download.profiles------------------------------------------------------------------------------------------------------------
+layer.name <- "wosis_latest_profiles"
+system.time(
+  ogr2ogr(src=wfs,
+          dst=wosis.dir.name,
+          layer=layer.name,
+          f = "ESRI Shapefile",
+          overwrite=TRUE,
+          skipfailures=TRUE)
+)
+
+
+## ----download.profiles.india------------------------------------------------------------------------------------------------------
+wosis.dir.name.india <- "./wosis_latest/india"
+if (!file.exists(wosis.dir.name.india)) dir.create(wosis.dir.name.india)
+layer.name <- "wosis_latest_profiles"
+system.time(
+  ogr2ogr(src=wfs,
+          dst=wosis.dir.name.india,
+          layer=layer.name,
+          f = "ESRI Shapefile",
+          where="country_name='India'",
+          overwrite=TRUE,
+          skipfailures=TRUE)
+)
+
+
+## ----import.profiles--------------------------------------------------------------------------------------------------------------
+layer.name <- "ms_wosis_latest_profiles"
+profiles <- readOGR(dsn=wosis.dir.name, layer=layer.name,
+                stringsAsFactors = FALSE)
+class(profiles)
+dim(profiles)
+names(profiles@data)
+head(profiles@data, 1)
+
+
+## ----import.profiles.india--------------------------------------------------------------------------------------------------------
+profiles.india <- readOGR(dsn=wosis.dir.name.india, layer=layer.name,
+                stringsAsFactors = FALSE)
+dim(profiles.india)
+
+
+## ----download.bd------------------------------------------------------------------------------------------------------------------
+layer.name <- "wosis_latest_bdfi33"
+system.time(
+  ogr2ogr(src=wfs,
+          dst=wosis.dir.name,
+          layer=layer.name,
+          f = "ESRI Shapefile",
+          overwrite=TRUE,
+          skipfailures=TRUE)
+)
+
+
+## ----download.bd.india------------------------------------------------------------------------------------------------------------
+system.time(
+  ogr2ogr(src=wfs,
+          dst=wosis.dir.name.india,
+          layer=layer.name,
+          f = "ESRI Shapefile",
+          where="country_name='India'",
+          overwrite=TRUE,
+          skipfailures=TRUE)
+)
+
+
+## ----import.bd--------------------------------------------------------------------------------------------------------------------
+# here strings are just that, not to be interpreted as R factors
+layer.name <- "ms_wosis_latest_bdfi33"
+bd33 <- readOGR(dsn=wosis.dir.name, layer=layer.name,
+                stringsAsFactors = FALSE)
+class(bd33)
+names(bd33@data)
+head(bd33@data, 1)
+
+
+## ----head.bd----------------------------------------------------------------------------------------------------------------------
+head(bd33@data$bdfi33_val)
+
+
+## ----example.bd-------------------------------------------------------------------------------------------------------------------
+bd33@data[75:80, c("profile_id","upper_dept","lower_dept","bdfi33_val","bdfi33_v_1")]
+
+
+## ----import.bd.india--------------------------------------------------------------------------------------------------------------
+layer.name <- "ms_wosis_latest_bdfi33"
+bd33.india <- readOGR(dsn=wosis.dir.name.india, layer=layer.name,
+                stringsAsFactors = FALSE)
+class(bd33.india)
+dim(bd33.india)
+(profile.id.india <- unique(bd33.india@data$profile_id))
+names(bd33.india@data)
+hist(bd33.india@data$bdfi33_v_1, main="Bulk density, Layers in India", 
+     xlab="g cm-3")
+
+
+## ---------------------------------------------------------------------------------------------------------------------------------
+
+
+
+## ----download.profiles.csv--------------------------------------------------------------------------------------------------------
+wosis.dir.name.ceu <- "./wosis_latest/central_europe"
+if (!file.exists(wosis.dir.name.ceu)) dir.create(wosis.dir.name.ceu)
+src.layer.name <- "ms:wosis_latest_profiles"
+dst.layer.name <- "wosis_latest_profiles_ceu"
+(dst.target.name <- paste0(wosis.dir.name.ceu,"/",dst.layer.name,".csv"))
+if (file.exists(dst.target.name)) { file.remove(dst.target.name) }
+ogr2ogr(src=wfs, 
+        dst=dst.target.name,
+        layer=src.layer.name,
+        f="CSV",
+        spat=c(6, 48, 8, 50),
+        overwrite=TRUE)
+round(file.info(dst.target.name)$size/1024,1)
+
+
+## ----import.profiles.csv----------------------------------------------------------------------------------------------------------
+layer.name <- "wosis_latest_profiles_ceu"
+system.time(
+  profiles.ceu <- read.csv(paste0(wosis.dir.name.ceu, "/",layer.name,".csv"),
+                stringsAsFactors = FALSE)
+)
+names(profiles.ceu)
+
+
+## ----download.bd.csv--------------------------------------------------------------------------------------------------------------
+src.layer.name <- "wosis_latest_bdfi33"
+dst.layer.name <- "wosis_latest_bdfi33"
+(dst.target.name <- paste0(wosis.dir.name,"/",dst.layer.name,".csv"))
+if (file.exists(dst.target.name)) { file.remove(dst.target.name) }
+ogr2ogr(src=wfs,
+        dst=dst.target.name,
+        layer=src.layer.name,
+        f="CSV",
+        overwrite=TRUE)
+file.info(dst.target.name)$size
+
+
+## ----import.bd.csv----------------------------------------------------------------------------------------------------------------
+layer.name <- "wosis_latest_bdfi33"
+system.time(
+  bd33.pts <- read.csv(file=paste0(wosis.dir.name,"/",layer.name,".csv"),
+                stringsAsFactors = FALSE)
+)
+dim(bd33.pts)
+names(bd33.pts)
+
+
+## ----download.profile.gpkg--------------------------------------------------------------------------------------------------------
+layer.name <- "wosis_latest_profiles"
+(dst.target.name <- paste0(wosis.dir.name,"/", layer.name, ".gpkg"))
+if (file.exists(dst.target.name)) { file.remove(dst.target.name) }
+system.time(
+  ogr2ogr(src=wfs,
+          dst=dst.target.name,
+          layer=layer.name,
+          f = "GPKG",
+          overwrite=TRUE,
+          skipfailures=TRUE)
+)
+
+
+## ----read.profile.gpkg------------------------------------------------------------------------------------------------------------
+profiles.gpkg <- st_read(dst.target.name)
+class(profiles.gpkg)
+dim(profiles.gpkg)
+names(profiles.gpkg)
+
+
+## ---------------------------------------------------------------------------------------------------------------------------------
+head(sort(table(profiles.gpkg$country_name), decreasing=TRUE))
+
+
+## ----map.profiles.gpkg, fig.width=12, fig.height=12-------------------------------------------------------------------------------
+ggplot(data=profiles.gpkg[(profiles.gpkg$country_name=="Brazil"), ]) +
+  aes(col=dataset_id) +
+  geom_sf(shape=21, size=0.8, fill="black")
+
+
+## ----coordinates.profiles---------------------------------------------------------------------------------------------------------
+coordinates(profiles.ceu) <- c("longitude", "latitude")
+proj4string(profiles.ceu) <- CRS("+init=epsg:4326")
+class(profiles.ceu)
+
+
+## ----coordinates.profiles.sf------------------------------------------------------------------------------------------------------
+profiles.ceu.sf <- st_as_sf(profiles.ceu)
+class(profiles.ceu.sf)
+
+
+## ----profiles.wrb-----------------------------------------------------------------------------------------------------------------
+table(profiles.ceu@data$cwrb_reference_soil_group)
+table(profiles.ceu.sf$cwrb_reference_soil_group)
+
+
+## ----map.profiles.wrb-------------------------------------------------------------------------------------------------------------
+ggplot(data=profiles.ceu.sf) +
+  aes(col=cwrb_reference_soil_group) +
+  geom_sf()
+
+
+## ----load.aqp---------------------------------------------------------------------------------------------------------------------
+library(aqp)            # Algorithms for Quantitative Pedology
+
+
+## ----bd.aqp-----------------------------------------------------------------------------------------------------------------------
+ds.aqp <- bd33@data
+depths(ds.aqp) <- profile_id ~ upper_dept + lower_dept
+is(ds.aqp)
+slotNames(ds.aqp)
+str(ds.aqp@site)
+str(ds.aqp@horizons)
+head(ds.aqp@site)
+head(ds.aqp@horizons[c(2,5,6,7,9)],12)
+
+
+## ----plot.bd.spc,fig.width=12, fig.height=8---------------------------------------------------------------------------------------
+plotSPC(ds.aqp[1:24,], name="layer_name", color='bdfi33_v_1')
+