Commit 693d03ae authored by Batjes, Niels's avatar Batjes, Niels
Browse files

Replace WoSIS_Latest_With_R.R

parent de768c1f
## ----setup, include=FALSE--------------------------------------------------------------------------------------------------------- ## ----setup, include=FALSE----------------------------------------------------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, warnings = FALSE, purl = FALSE) knitr::opts_chunk$set(echo = TRUE, warnings = FALSE, purl = FALSE)
## ----load.package----------------------------------------------------------------------------------------------------------------- ## ----load.package------------------------------------------------------------------------------------------------------
library(rgdal) # GDAL access from R library(rgdal) # GDAL access from R
library(gdalUtils) # wrappers for GDAL utility programs that could be library(gdalUtils) # wrappers for GDAL utility programs that could be
# called from the command line # called from the command line
...@@ -13,7 +13,7 @@ library(ggplot2) # gpplot graphics ...@@ -13,7 +13,7 @@ library(ggplot2) # gpplot graphics
## --------------------------------------------------------------------------------------------------------------------------------- ## ----------------------------------------------------------------------------------------------------------------------
drivers <- ogrDrivers() drivers <- ogrDrivers()
# print(drivers) # print(drivers)
ix <- grep("GPKG", drivers$name, fixed=TRUE) ix <- grep("GPKG", drivers$name, fixed=TRUE)
...@@ -22,18 +22,18 @@ ix <- grep("ESRI", drivers$name, fixed=TRUE) ...@@ -22,18 +22,18 @@ ix <- grep("ESRI", drivers$name, fixed=TRUE)
drivers[ix,] drivers[ix,]
## ----specify.wfs.address---------------------------------------------------------------------------------------------------------- ## ----specify.wfs.address-----------------------------------------------------------------------------------------------
wfs <- "WFS:http://maps.isric.org/mapserv?map=/map/wosis_latest.map" wfs <- "WFS:https://maps.isric.org/mapserv?map=/map/wosis_latest.map"
## ----ogrinfo.wfs------------------------------------------------------------------------------------------------------------------ ## ----ogrinfo.wfs-------------------------------------------------------------------------------------------------------
system.time( system.time(
layers.info <- ogrinfo(wfs, ro=TRUE, so=TRUE, q=TRUE) layers.info <- ogrinfo(wfs, ro=TRUE, so=TRUE, q=TRUE)
) )
cat(layers.info, sep = "\n") cat(layers.info, sep = "\n")
## ----ogrinfo.site----------------------------------------------------------------------------------------------------------------- ## ----ogrinfo.site------------------------------------------------------------------------------------------------------
system.time( system.time(
profiles.info <- profiles.info <-
ogrinfo(wfs, layer="ms:wosis_latest_profiles", ogrinfo(wfs, layer="ms:wosis_latest_profiles",
...@@ -42,7 +42,7 @@ system.time( ...@@ -42,7 +42,7 @@ system.time(
cat(profiles.info, sep="\n") cat(profiles.info, sep="\n")
## ----ogrinfo.eu.profiles.1-------------------------------------------------------------------------------------------------------- ## ----ogrinfo.eu.profiles.1---------------------------------------------------------------------------------------------
system.time( system.time(
central.eu.profiles.info <- central.eu.profiles.info <-
ogrinfo(wfs, ro=TRUE, so=TRUE, q=FALSE, ogrinfo(wfs, ro=TRUE, so=TRUE, q=FALSE,
...@@ -53,7 +53,7 @@ head(central.eu.profiles.info, 8) ...@@ -53,7 +53,7 @@ head(central.eu.profiles.info, 8)
central.eu.profiles.info[41:60] central.eu.profiles.info[41:60]
## ----ogrinfo.eu.profiles.2-------------------------------------------------------------------------------------------------------- ## ----ogrinfo.eu.profiles.2---------------------------------------------------------------------------------------------
ix.f <- grep("Feature Count", central.eu.profiles.info) ix.f <- grep("Feature Count", central.eu.profiles.info)
central.eu.profiles.info[ix.f] central.eu.profiles.info[ix.f]
ix.e <- grep("Extent", central.eu.profiles.info) ix.e <- grep("Extent", central.eu.profiles.info)
...@@ -62,13 +62,13 @@ ix.g <- grep("GEOGCRS", central.eu.profiles.info) ...@@ -62,13 +62,13 @@ ix.g <- grep("GEOGCRS", central.eu.profiles.info)
cat(paste(central.eu.profiles.info[ix.g:(ix.g+17)], collapse="\n")) cat(paste(central.eu.profiles.info[ix.g:(ix.g+17)], collapse="\n"))
## ----show.fields------------------------------------------------------------------------------------------------------------------ ## ----show.fields-------------------------------------------------------------------------------------------------------
ix.p <- grep("Geometry Column", central.eu.profiles.info) ix.p <- grep("Geometry Column", central.eu.profiles.info)
n <- length(central.eu.profiles.info) n <- length(central.eu.profiles.info)
central.eu.profiles.info[ix.p+1:n] central.eu.profiles.info[ix.p+1:n]
## ----ogrinfo.india.profiles------------------------------------------------------------------------------------------------------- ## ----ogrinfo.india.profiles--------------------------------------------------------------------------------------------
system.time( system.time(
india.profiles.info <- india.profiles.info <-
ogrinfo(wfs, ro=TRUE, so=TRUE, q=FALSE, ogrinfo(wfs, ro=TRUE, so=TRUE, q=FALSE,
...@@ -77,14 +77,14 @@ system.time( ...@@ -77,14 +77,14 @@ system.time(
) )
## ----ogrinfo.india.profiles.2----------------------------------------------------------------------------------------------------- ## ----ogrinfo.india.profiles.2------------------------------------------------------------------------------------------
ix.f <- grep("Feature Count", india.profiles.info) ix.f <- grep("Feature Count", india.profiles.info)
india.profiles.info[ix.f] india.profiles.info[ix.f]
ix.e <- grep("Extent", india.profiles.info) ix.e <- grep("Extent", india.profiles.info)
india.profiles.info[ix.e] india.profiles.info[ix.e]
## ----ogrinfo.properties, warnings------------------------------------------------------------------------------------------------- ## ----ogrinfo.properties, warnings--------------------------------------------------------------------------------------
system.time( system.time(
property.info <- ogrinfo(wfs, layer="ms:wosis_latest_bdfi33", property.info <- ogrinfo(wfs, layer="ms:wosis_latest_bdfi33",
ro=TRUE, so=TRUE, q=FALSE) ro=TRUE, so=TRUE, q=FALSE)
...@@ -94,7 +94,7 @@ ix.f <- grep("Feature Count", property.info) ...@@ -94,7 +94,7 @@ ix.f <- grep("Feature Count", property.info)
property.info[ix.f] property.info[ix.f]
## ----ogrinfo.bd.india------------------------------------------------------------------------------------------------------------- ## ----ogrinfo.bd.india--------------------------------------------------------------------------------------------------
system.time( system.time(
bd.india.info <- ogrinfo(wfs, layer="ms:wosis_latest_bdfi33", bd.india.info <- ogrinfo(wfs, layer="ms:wosis_latest_bdfi33",
where="country_name='India' AND upper_depth > 100", where="country_name='India' AND upper_depth > 100",
...@@ -106,12 +106,12 @@ bd.india.info[ix.f] ...@@ -106,12 +106,12 @@ bd.india.info[ix.f]
split=": ", fixed=TRUE)[[1]][2])) split=": ", fixed=TRUE)[[1]][2]))
## ----setup.local.dir-------------------------------------------------------------------------------------------------------------- ## ----setup.local.dir---------------------------------------------------------------------------------------------------
wosis.dir.name <- "./wosis_latest" wosis.dir.name <- "./wosis_latest"
if (!file.exists(wosis.dir.name)) dir.create(wosis.dir.name) if (!file.exists(wosis.dir.name)) dir.create(wosis.dir.name)
## ----download.profiles------------------------------------------------------------------------------------------------------------ ## ----download.profiles-------------------------------------------------------------------------------------------------
layer.name <- "wosis_latest_profiles" layer.name <- "wosis_latest_profiles"
system.time( system.time(
ogr2ogr(src=wfs, ogr2ogr(src=wfs,
...@@ -123,7 +123,7 @@ system.time( ...@@ -123,7 +123,7 @@ system.time(
) )
## ----download.profiles.india------------------------------------------------------------------------------------------------------ ## ----download.profiles.india-------------------------------------------------------------------------------------------
wosis.dir.name.india <- "./wosis_latest/india" wosis.dir.name.india <- "./wosis_latest/india"
if (!file.exists(wosis.dir.name.india)) dir.create(wosis.dir.name.india) if (!file.exists(wosis.dir.name.india)) dir.create(wosis.dir.name.india)
layer.name <- "wosis_latest_profiles" layer.name <- "wosis_latest_profiles"
...@@ -138,7 +138,7 @@ system.time( ...@@ -138,7 +138,7 @@ system.time(
) )
## ----import.profiles-------------------------------------------------------------------------------------------------------------- ## ----import.profiles---------------------------------------------------------------------------------------------------
layer.name <- "ms_wosis_latest_profiles" layer.name <- "ms_wosis_latest_profiles"
profiles <- readOGR(dsn=wosis.dir.name, layer=layer.name, profiles <- readOGR(dsn=wosis.dir.name, layer=layer.name,
stringsAsFactors = FALSE) stringsAsFactors = FALSE)
...@@ -148,13 +148,13 @@ names(profiles@data) ...@@ -148,13 +148,13 @@ names(profiles@data)
head(profiles@data, 1) head(profiles@data, 1)
## ----import.profiles.india-------------------------------------------------------------------------------------------------------- ## ----import.profiles.india---------------------------------------------------------------------------------------------
profiles.india <- readOGR(dsn=wosis.dir.name.india, layer=layer.name, profiles.india <- readOGR(dsn=wosis.dir.name.india, layer=layer.name,
stringsAsFactors = FALSE) stringsAsFactors = FALSE)
dim(profiles.india) dim(profiles.india)
## ----download.bd------------------------------------------------------------------------------------------------------------------ ## ----download.bd-------------------------------------------------------------------------------------------------------
layer.name <- "wosis_latest_bdfi33" layer.name <- "wosis_latest_bdfi33"
system.time( system.time(
ogr2ogr(src=wfs, ogr2ogr(src=wfs,
...@@ -166,7 +166,7 @@ system.time( ...@@ -166,7 +166,7 @@ system.time(
) )
## ----download.bd.india------------------------------------------------------------------------------------------------------------ ## ----download.bd.india-------------------------------------------------------------------------------------------------
system.time( system.time(
ogr2ogr(src=wfs, ogr2ogr(src=wfs,
dst=wosis.dir.name.india, dst=wosis.dir.name.india,
...@@ -178,7 +178,7 @@ system.time( ...@@ -178,7 +178,7 @@ system.time(
) )
## ----import.bd-------------------------------------------------------------------------------------------------------------------- ## ----import.bd---------------------------------------------------------------------------------------------------------
# here strings are just that, not to be interpreted as R factors # here strings are just that, not to be interpreted as R factors
layer.name <- "ms_wosis_latest_bdfi33" layer.name <- "ms_wosis_latest_bdfi33"
bd33 <- readOGR(dsn=wosis.dir.name, layer=layer.name, bd33 <- readOGR(dsn=wosis.dir.name, layer=layer.name,
...@@ -188,15 +188,15 @@ names(bd33@data) ...@@ -188,15 +188,15 @@ names(bd33@data)
head(bd33@data, 1) head(bd33@data, 1)
## ----head.bd---------------------------------------------------------------------------------------------------------------------- ## ----head.bd-----------------------------------------------------------------------------------------------------------
head(bd33@data$bdfi33_val) head(bd33@data$bdfi33_val)
## ----example.bd------------------------------------------------------------------------------------------------------------------- ## ----example.bd--------------------------------------------------------------------------------------------------------
bd33@data[75:80, c("profile_id","upper_dept","lower_dept","bdfi33_val","bdfi33_v_1")] bd33@data[75:80, c("profile_id","upper_dept","lower_dept","bdfi33_val","bdfi33_v_1")]
## ----import.bd.india-------------------------------------------------------------------------------------------------------------- ## ----import.bd.india---------------------------------------------------------------------------------------------------
layer.name <- "ms_wosis_latest_bdfi33" layer.name <- "ms_wosis_latest_bdfi33"
bd33.india <- readOGR(dsn=wosis.dir.name.india, layer=layer.name, bd33.india <- readOGR(dsn=wosis.dir.name.india, layer=layer.name,
stringsAsFactors = FALSE) stringsAsFactors = FALSE)
...@@ -208,11 +208,7 @@ hist(bd33.india@data$bdfi33_v_1, main="Bulk density, Layers in India", ...@@ -208,11 +208,7 @@ hist(bd33.india@data$bdfi33_v_1, main="Bulk density, Layers in India",
xlab="g cm-3") xlab="g cm-3")
## --------------------------------------------------------------------------------------------------------------------------------- ## ----download.profiles.csv---------------------------------------------------------------------------------------------
## ----download.profiles.csv--------------------------------------------------------------------------------------------------------
wosis.dir.name.ceu <- "./wosis_latest/central_europe" wosis.dir.name.ceu <- "./wosis_latest/central_europe"
if (!file.exists(wosis.dir.name.ceu)) dir.create(wosis.dir.name.ceu) if (!file.exists(wosis.dir.name.ceu)) dir.create(wosis.dir.name.ceu)
src.layer.name <- "ms:wosis_latest_profiles" src.layer.name <- "ms:wosis_latest_profiles"
...@@ -228,7 +224,7 @@ ogr2ogr(src=wfs, ...@@ -228,7 +224,7 @@ ogr2ogr(src=wfs,
round(file.info(dst.target.name)$size/1024,1) round(file.info(dst.target.name)$size/1024,1)
## ----import.profiles.csv---------------------------------------------------------------------------------------------------------- ## ----import.profiles.csv-----------------------------------------------------------------------------------------------
layer.name <- "wosis_latest_profiles_ceu" layer.name <- "wosis_latest_profiles_ceu"
system.time( system.time(
profiles.ceu <- read.csv(paste0(wosis.dir.name.ceu, "/",layer.name,".csv"), profiles.ceu <- read.csv(paste0(wosis.dir.name.ceu, "/",layer.name,".csv"),
...@@ -237,7 +233,7 @@ system.time( ...@@ -237,7 +233,7 @@ system.time(
names(profiles.ceu) names(profiles.ceu)
## ----download.bd.csv-------------------------------------------------------------------------------------------------------------- ## ----download.bd.csv---------------------------------------------------------------------------------------------------
src.layer.name <- "wosis_latest_bdfi33" src.layer.name <- "wosis_latest_bdfi33"
dst.layer.name <- "wosis_latest_bdfi33" dst.layer.name <- "wosis_latest_bdfi33"
(dst.target.name <- paste0(wosis.dir.name,"/",dst.layer.name,".csv")) (dst.target.name <- paste0(wosis.dir.name,"/",dst.layer.name,".csv"))
...@@ -247,10 +243,10 @@ ogr2ogr(src=wfs, ...@@ -247,10 +243,10 @@ ogr2ogr(src=wfs,
layer=src.layer.name, layer=src.layer.name,
f="CSV", f="CSV",
overwrite=TRUE) overwrite=TRUE)
file.info(dst.target.name)$size file.info(dst.target.name)$size/1024/1024
## ----import.bd.csv---------------------------------------------------------------------------------------------------------------- ## ----import.bd.csv-----------------------------------------------------------------------------------------------------
layer.name <- "wosis_latest_bdfi33" layer.name <- "wosis_latest_bdfi33"
system.time( system.time(
bd33.pts <- read.csv(file=paste0(wosis.dir.name,"/",layer.name,".csv"), bd33.pts <- read.csv(file=paste0(wosis.dir.name,"/",layer.name,".csv"),
...@@ -260,7 +256,7 @@ dim(bd33.pts) ...@@ -260,7 +256,7 @@ dim(bd33.pts)
names(bd33.pts) names(bd33.pts)
## ----download.profile.gpkg-------------------------------------------------------------------------------------------------------- ## ----download.profile.gpkg---------------------------------------------------------------------------------------------
layer.name <- "wosis_latest_profiles" layer.name <- "wosis_latest_profiles"
(dst.target.name <- paste0(wosis.dir.name,"/", layer.name, ".gpkg")) (dst.target.name <- paste0(wosis.dir.name,"/", layer.name, ".gpkg"))
if (file.exists(dst.target.name)) { file.remove(dst.target.name) } if (file.exists(dst.target.name)) { file.remove(dst.target.name) }
...@@ -272,52 +268,53 @@ system.time( ...@@ -272,52 +268,53 @@ system.time(
overwrite=TRUE, overwrite=TRUE,
skipfailures=TRUE) skipfailures=TRUE)
) )
file.info(dst.target.name)$size/1024/1024
## ----read.profile.gpkg------------------------------------------------------------------------------------------------------------ ## ----read.profile.gpkg-------------------------------------------------------------------------------------------------
profiles.gpkg <- st_read(dst.target.name) profiles.gpkg <- st_read(dst.target.name)
class(profiles.gpkg) class(profiles.gpkg)
dim(profiles.gpkg) dim(profiles.gpkg)
names(profiles.gpkg) names(profiles.gpkg)
## --------------------------------------------------------------------------------------------------------------------------------- ## ----------------------------------------------------------------------------------------------------------------------
head(sort(table(profiles.gpkg$country_name), decreasing=TRUE)) head(sort(table(profiles.gpkg$country_name), decreasing=TRUE))
## ----map.profiles.gpkg, fig.width=12, fig.height=12------------------------------------------------------------------------------- ## ----map.profiles.gpkg, fig.width=12, fig.height=12--------------------------------------------------------------------
ggplot(data=profiles.gpkg[(profiles.gpkg$country_name=="Brazil"), ]) + ggplot(data=profiles.gpkg[(profiles.gpkg$country_name=="Brazil"), ]) +
aes(col=dataset_id) + aes(col=dataset_id) +
geom_sf(shape=21, size=0.8, fill="black") geom_sf(shape=21, size=0.8, fill="black")
## ----coordinates.profiles--------------------------------------------------------------------------------------------------------- ## ----coordinates.profiles----------------------------------------------------------------------------------------------
coordinates(profiles.ceu) <- c("longitude", "latitude") coordinates(profiles.ceu) <- c("longitude", "latitude")
proj4string(profiles.ceu) <- CRS("+init=epsg:4326") proj4string(profiles.ceu) <- CRS("+init=epsg:4326")
class(profiles.ceu) class(profiles.ceu)
## ----coordinates.profiles.sf------------------------------------------------------------------------------------------------------ ## ----coordinates.profiles.sf-------------------------------------------------------------------------------------------
profiles.ceu.sf <- st_as_sf(profiles.ceu) profiles.ceu.sf <- st_as_sf(profiles.ceu)
class(profiles.ceu.sf) class(profiles.ceu.sf)
## ----profiles.wrb----------------------------------------------------------------------------------------------------------------- ## ----profiles.wrb------------------------------------------------------------------------------------------------------
table(profiles.ceu@data$cwrb_reference_soil_group) table(profiles.ceu@data$cwrb_reference_soil_group)
table(profiles.ceu.sf$cwrb_reference_soil_group) table(profiles.ceu.sf$cwrb_reference_soil_group)
## ----map.profiles.wrb------------------------------------------------------------------------------------------------------------- ## ----map.profiles.wrb--------------------------------------------------------------------------------------------------
ggplot(data=profiles.ceu.sf) + ggplot(data=profiles.ceu.sf) +
aes(col=cwrb_reference_soil_group) + aes(col=cwrb_reference_soil_group) +
geom_sf() geom_sf()
## ----load.aqp--------------------------------------------------------------------------------------------------------------------- ## ----load.aqp----------------------------------------------------------------------------------------------------------
library(aqp) # Algorithms for Quantitative Pedology library(aqp) # Algorithms for Quantitative Pedology
## ----bd.aqp----------------------------------------------------------------------------------------------------------------------- ## ----bd.aqp------------------------------------------------------------------------------------------------------------
ds.aqp <- bd33@data ds.aqp <- bd33@data
depths(ds.aqp) <- profile_id ~ upper_dept + lower_dept depths(ds.aqp) <- profile_id ~ upper_dept + lower_dept
is(ds.aqp) is(ds.aqp)
...@@ -328,6 +325,6 @@ head(ds.aqp@site) ...@@ -328,6 +325,6 @@ head(ds.aqp@site)
head(ds.aqp@horizons[c(2,5,6,7,9)],12) head(ds.aqp@horizons[c(2,5,6,7,9)],12)
## ----plot.bd.spc,fig.width=12, fig.height=8--------------------------------------------------------------------------------------- ## ----plot.bd.spc,fig.width=12, fig.height=8----------------------------------------------------------------------------
plotSPC(ds.aqp[1:24,], name="layer_name", color='bdfi33_v_1') plotSPC(ds.aqp[1:24,], name="layer_name", color='bdfi33_v_1')
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment