André A. Stella†, João P. S. Pavan†, Mauricio S. Araújo¹, D Bruno F. Fregonezi,, Innocent Vulou Unzimai, Erica P. Leles, Michelle F. Santos, Peter Goldsmith, Godfree Chigeza, Brian W. Diers, Tony Gathungu, ÂJames Njoroge, Maria I. Zucchi³˒⁴, José B. Pinheiro † These authors contributed equally to this work.
Stella, A.A., Pavan, J.P.S., Araújo, M.S. et al. Soybean selection in Kenya enhanced by multi-trait and genotype-by-environment interaction modeling. Sci Rep 15, 27575 (2025). https://doi.org/10.1038/s41598-025-10654-2
library(tmap)
library(dplyr)
library(rnaturalearth)
library(ggplot2)
library(sf)
library(terra)
library(rnaturalearthdata)
library(ggspatial)
library(patchwork)
library(dplyr)
library(forcats)
library(viridis)
library(RColorBrewer)
library(nasapower)
library(raster)
library(terra)
library(SoilType)
library(soiltexture)
library(tidyterra)
library(geodata)
library(raster)
library(sf)
library(ggplot2)
library(rnaturalearth)
library(rnaturalearthdata)
library(cowplot)
library(viridis)
library(elevatr)
library(ggspatial)
library(ggpubr)
library(here)raster_agro33 = rast(here("Raster_Agro33", "agro33.tif"))
#Choose Africa
afr = ne_countries(scale = "medium", continent = "Africa", returnclass = "sf")africabbox= st_bbox(afr)
agro_afr = crop(raster_agro33, ext(africabbox))agrodf = as.data.frame(agro_afr, xy= T)
colnames(agrodf) = c("lon", "lat", "AEZ_class")
kenya = ne_countries(scale = "medium", country = "Kenya", returnclass = "sf")kenya_env = readxl::read_xlsx(here("Kenya_regions", "regkenya.xlsx"))afrshp = shapefile(here("Africa_Biomes_Dataset_shapefile", "Africa_Biomes_Dataset.shp"))
afrshp_sf = sf::st_as_sf(afrshp)
kenya_sf= st_as_sf(kenya)
if (st_crs(afrshp_sf) != st_crs(kenya_sf)) {
kenya_sf = st_transform(kenya_sf, st_crs(afrshp_sf))
}
kenyabio = st_intersection(afrshp_sf, kenya_sf)themekenya = theme_minimal(base_size = 14) +
theme(
legend.position = "right",
legend.title = element_text(size = 12),
legend.text = element_text(size = 10),
axis.title = element_text(size = 12),
axis.text = element_text(size = 12),
plot.margin = margin(0,0,0,0, "cm")
)biomesfrica = ggplot() +
geom_sf(data = afrshp_sf, aes(fill = ECO_NAME),
color = "black", alpha = 0.6) +
geom_sf(data = kenya, fill = NA, color = "red", size = 6) +
scale_fill_manual(name = "Biomes",
values = c("#477036", "#82b74b","#f0ead6",
"#96ded1", "#adcca8", "#996515", "#89922e",
"#a87c6d", "#aa98a9", "#bc5a45","#e3eaa7")) +
theme(legend.position = "left") +
annotation_scale(location = "bl", width_hint = 0.3) +
annotation_north_arrow(location = "tr", which_north = "true",
style = north_arrow_fancy_orienteering,
height = unit(0.8, "cm"), width = unit(0.8, "cm")) +
labs(x= "Longitude", y = "Latitude") +
theme_minimal(base_size = 14) +
themekenya
Biomes_Africa = biomesfrica + ggtitle("A")
biomesfrica
kenya_env = kenya_env %>%
mutate(lat = as.numeric(lat),
long = as.numeric(long))
kenya_envsf = st_as_sf(kenya_env, coords = c("long", "lat"), crs = 4326)
kenya_envsf <- st_transform(kenya_envsf, crs = 4326)
env_biomeskenya = ggplot() +
geom_sf(data = kenyabio, aes(fill=ECO_NAME), color = "black", alpha = 1) +
geom_sf(data = kenya, fill = NA, color = "black", size = 1.2) +
scale_fill_manual(name = "Biomes",
values = c("#477036","#82b74b", "#f0ead6",
"#96ded1","#996515","#aa98a9","#e3eaa7"),
guide = guide_legend(order = 1)) +
geom_sf(data = kenya_envsf, aes(shape = env), size = 4, color = "black") +
scale_shape_manual(name = "Region", values = c(15:18, 8:12)) +
themekenya +
annotation_scale(location = "bl", width_hint = 0.3) +
annotation_north_arrow(location = "tr", which_north = "true",
style = north_arrow_fancy_orienteering,
height = unit(0.8, "cm"), width = unit(0.8, "cm")) +
labs(x= "Longitude", y = "Latitude")
env_biomeskenya
Biomes_Kenya = env_biomeskenya + ggtitle("B")soilraster = rast(here("Soil_Raster", "hwsd_domi.tif"))
kenya_soilraster = crop(soilraster, vect(kenya))
kenya_soilraster= mask(kenya_soilraster, vect(kenya))
soilclasses = unique(values(kenya_soilraster))
soilclasses = na.omit(soilclasses)soilcolors = c("#ed872d", "#be6516", "#fbc97f", "#f5f3ce", "#efd25c",
"#b2beb5", "#4f372d", "#928f88", "#d0674f", "#8d4338",
"#eae679", "#c28a44", "#f99379", "#f3e5ab", "#e8dad4",
"#e3c882", "#acadad", "#d7c0b6","#e96957", "#ac826d",
"#c8c8a7", "#aa381e","#f5f3ce","#deba92", "#727a77")soil_class_names <- data.frame(
soil_class = c(7, 9, 27, 26, 13, 4, 22, 25, 8, 31, 16, 18, 17,
10, 28, 21, 2, 11, 1, 19, 14, 6, 3, 12, 32),
class_name = c(
"Latosol", "Argisol", "Cambisol", "Chernozem", "Ferralsol",
"Gleysol", "Histosol", "Leptosol", "Luvisol", "Nitisol",
"Phaeozem", "Planosol", "Podzol", "Regosol", "Solonchak",
"Solonetz", "Vertisol", "Arenosol", "Lixisol", "Kastanozem",
"Calcisol", "Gypsisol", "Durisol", "Plinthosol", "Fluvisol"
)
)
soildf = as.data.frame(kenya_soilraster, xy = T)
colnames(soildf) = c("x", "y", "soil_class")
soildf = na.omit(soildf)
soildf= merge(soildf, soil_class_names, by = "soil_class")soilmapkenya =
ggplot() +
geom_raster(data = soildf, aes(x = x, y = y, fill = class_name)) +
geom_sf(data = kenya, fill = NA, color = "black", size = 1.2) +
geom_sf(data = kenya_envsf, aes(shape = env), size = 4, color = "black") +
scale_fill_manual(name = "Soil Classes", values = soilcolors) +
scale_shape_manual(name = "Region", values = c(15:18, 8:12)) +
coord_sf(xlim = c(34, 42), ylim = c(-5, 6), expand = TRUE) +
themekenya +
labs(
x = "Longitude",
y = "Latitude") +
annotation_scale(location = "bl", width_hint = 0.2,
bar_cols = c("black", "white"), text_cex = 0.8,
unit_category = "metric",
line_width = 0.5, text_col = "black") +
annotation_north_arrow(
location = "tr",
which_north = "true",
style = north_arrow_fancy_orienteering(),
height = unit(0.8, "cm"),
width = unit(0.8, "cm")
)
soilmapkenya
Soil_Kenya = soilmapkenya + ggtitle("C")kenya = ne_countries(scale = "medium", country = "Kenya", returnclass = "sf")30 s resolution (~1km)
elev_kenya = elevation_30s(country = "Kenya", path = tempdir())
#elev_kenya = na.omit(elev_kenya)ele_kenyadf = as.data.frame(elev_kenya, xy = TRUE)colnames(ele_kenyadf)[3] = "altitude"kenya = ne_countries(scale = "medium", country = "Kenya", returnclass = "sf")
kenyabox = st_bbox(kenya)
longseq = seq(kenyabox$xmin, kenyabox$xmax, length.out = 20)
latseq = seq(kenyabox$ymin, kenyabox$ymax, length.out=20)
gridpoints = expand.grid(longseq, latseq)
colnames(gridpoints) = c("Long", "Lat")
kenyapoints = st_as_sf(gridpoints, coords = c("Long", "Lat"), crs = 4326)
altitudekenya = get_elev_raster(
locations = kenyapoints, z=7, clip = "bbox", src = "aws")
altituderast = rast(altitudekenya)
kenyavect = vect(kenya)
altitudekenyacrop = crop(altituderast, kenyavect)
altitudekenyamasked = mask(altitudekenyacrop, kenyavect)
altitudekenyadf = as.data.frame(altitudekenyamasked, xy = T)
colnames(altitudekenyadf) = c("Long", "Lat", "Elevation")cp1 = c("#f0e6d2","#edf08e","#66bb6a","#1b5e20", "#827717",
"#8d6e63", "#795548", "#6d4c41","#4e342e", "#7e57c2","#dcd0ff")
common_limits = coord_sf(xlim = c(33,42), ylim = c(-5,6), expand = F)kenyaaltitude = ggplot() +
geom_tile(data = altitudekenyadf, aes(x = Long, y = Lat, fill = Elevation)) +
geom_sf(data = kenya, fill = NA, color = "black", size = 0.5) +
scale_fill_gradientn(colors = cp1, name = "Altitude (m)") +
annotation_scale(location = "bl", width_hint = 0.3) +
annotation_north_arrow(location = "tr", which_north = "true",
style = north_arrow_fancy_orienteering,
height = unit(0.8, "cm"), width = unit(0.8, "cm")) +
labs(x= "Longitude", y = "Latitude") +
geom_sf(data = kenya_envsf, aes(shape = env), size = 4, color = "black") +
scale_shape_manual(name = "Region", values = c(15:18, 8:12)) +
themekenya + common_limits
kenyaaltitude
Kenya_altitude = kenyaaltitude + ggtitle("D")mkenyalegend = ggplot() +
geom_sf(data = kenya_envsf, aes(shape = env), size = 4 , color = "black") +
scale_shape_manual(name = "Location", values = c (15:18, 8:12),
guide = guide_legend(nrow = 2)) +
theme_void() +
theme(legend.position = "bottom", legend.direction = "horizontal",
legend.title = element_text(size = 12),
legend.text = element_text(size = 10),
legend.margin = margin(t = 0, b = 0),
plot.margin = margin(0, 0, 0, 0)) ####Remove Location legend of the maps below
env_biomeskenya_noreg = env_biomeskenya + guides(shape = "none") +ggtitle("C")
soilmapkenya_noreg = soilmapkenya + guides(shape = "none") + ggtitle("B")
kenyaaltitude_noreg = kenyaaltitude + guides(shape ="none") + ggtitle("D")library(cowplot)
reglegendonly = get_legend(mkenyalegend)
kmap = (Biomes_Africa |soilmapkenya_noreg) /
(env_biomeskenya_noreg| kenyaaltitude_noreg)
kmap2= kmap / wrap_elements(reglegendonly)
kmapScale on map varies by more than 10%, scale bar may be inaccurate
Save
ggsave("kenyafinalmap2.pdf", width = 14, height = 10, dpi = 600)Have questions, want to collaborate, or found a bug?
Feel free to contact:
**João P. S. Pavan ** joaosilvapavan@usp.br