How to plot two raster results in the same graph with the legend of non-numeric variables？
I plot two graphs with the code below:
library(sf)
library(raster)
library(tidyverse)
library(RColorBrewer)
library(gridExtra)
library(grid)
plot_simulation <- function(input_file) {
Results <- read.csv(file = input_file, sep = ",", dec = ".", header = TRUE)
Results$CLocation <- gsub(" ", "", Results$wlieu) # 去除列 wlieu 的空格字符
medianYield <- aggregate(mafruit ~ CLocation, data = Results, FUN = median)
SimulatedPoints <- read.csv(file = "./LinPoints.csv", sep = ",", dec = ".", header = TRUE)
SimulatedPoints$CLocation <- gsub("^(\\d{2})_(\\d{1,2})$", "\\1_0\\2_v3test", SimulatedPoints$CLocation)
SimulatedPoints$CLocation <- gsub("^(\\d{2})_(\\d{3})$", "\\1_\\2_v3test", SimulatedPoints$CLocation)
SimulatedPointsCoords <- SimulatedPoints[, c("Grid_Code", "CLocation", "latitude", "longitude")]
SimulatedPointsCoords$CLocation <- gsub(" ", "", SimulatedPointsCoords$CLocation)
Results <- merge(SimulatedPointsCoords, medianYield, by = "CLocation", all = TRUE)
sf_Results <- st_as_sf(x = Results, coords = c("longitude", "latitude"), crs = 4326)
cell_size <- 30000
utm_crs <- st_crs(st_transform(sf_Results, 32631))
sf_Results_projected <- st_transform(sf_Results, utm_crs)
bounds <- st_bbox(sf_Results_projected)
x_min <- bounds["xmin"] - cell_size / 2
x_max <- bounds["xmax"] + cell_size / 2
y_min <- bounds["ymin"] - cell_size / 2
y_max <- bounds["ymax"] + cell_size / 2
raster_template <- raster(
xmn = x_min, xmx = x_max, ymn = y_min, ymx = y_max,
res = cell_size, crs = utm_crs$wkt
)
raster_Results <- rasterize(
x = st_coordinates(sf_Results_projected),
y = raster_template,
field = sf_Results$mafruit,
fun = median
)
url <- "https://raw.githubusercontent.com/gregoiredavid/france-geojson/master/regions.geojson"
france <- st_read(url)
france <- st_transform(france, crs = utm_crs) # 转换为 UTM 坐标系
raster_Results_wgs84 <- projectRaster(raster_Results, crs = "+proj=longlat +datum=WGS84")
france_wgs84 <- st_transform(france, crs = 4326)
color_breaks <- c(0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0)
#colors <- c("#E31A1C", "#FF7F00", "#FDBF6F", "#E9E4A6", "#A4D4A9", "#B2DF8A", "#33A02C", "#1F78B4")
colors <- c("#9B1A1A","#E31A1C", "#FF7F00", "#FDBF6F", "#B2DF8A", "#33A02C", "#1F78B4","#084E6A")
p <- ggplot() +
geom_sf(data = france_wgs84, fill = "white", color = "black", linewidth = 0.7) +
geom_tile(
data = as.data.frame(raster_Results_wgs84, xy = TRUE),
aes(x = x, y = y, fill = layer)
) +
scale_fill_gradientn(
colors = colors,
breaks = color_breaks,
labels = as.character(color_breaks),
name = expression("Yield (t·ha"^-1*")"),
na.value = NA,
limits = c(min(color_breaks), max(color_breaks))
) +
coord_sf(
xlim = c(st_bbox(france_wgs84)$xmin, st_bbox(france_wgs84)$xmax),
ylim = c(st_bbox(france_wgs84)$ymin, st_bbox(france_wgs84)$ymax),
expand = FALSE
) +
theme_minimal() +
theme(
plot.title = element_text(size = 12, face = "bold"),
legend.title = element_text(size = 10),
legend.text = element_text(size = 8),
legend.key.size = unit(1, "cm"),
legend.spacing.x = unit(0.5, "cm"),
legend.spacing.y = unit(0.5, "cm"),
panel.background = element_rect(fill = "white", color = "white"),
plot.background = element_rect(fill = "white", color = "white")
) +
guides(fill = guide_colorbar(
title = expression("Yield\n(t·ha"^-1*")"),
barwidth = 0.5,
barheight = 10,
title.position = "top",
title.hjust = 0.5,
label.position = "right",
direction = "vertical",
frame.colour = "black"
))
return(p)
}
files <- c(
"MGI_Pallador_irri_0.65_late_sowing.csv",
"MGI_Pallador_no_irri_0.65_late_sowing.csv"
)
plots <- list()
for (i in 1:length(files)) {
plots[[i]] <- plot_simulation(files[i])
}
grid.arrange(grobs = plots, ncol = 2, nrow = 1)
And I want the graph below:
But I want to combine these two figures into one figure:
Step 1.Combine the values in the same cell and choose the highest one.
Step 2.Assign two colors to these two files, and replace the highest value selected in Step 1 with the color corresponding to the file it came from. Additionally, create a legend using these two files.
Here is the expected result image:
"
"
IMPORTANT: The result map also needs to be raster type!
Original data:
LinPoints <- read.csv("LinPoints.csv")
dput(LinPoints)
structure(list(Grid_Code = c(31351652L, 31331654L, 31331655L,
31331657L, 31331658L, 31291905L, 31261906L, 31121911L), CLocation = c("73_63",
"73_63", "73_63", "73_63", "73_63", "73_73", "73_73", "73_73"
), elevation = c(1L, 1L, 1L, 1L, 1L, 20L, 20L, 20L), latitude = c(43.3697,
43.3909, 43.3926, 43.3961, 43.3978, 43.8155, 43.8436, 43.9751
), longitude = c(-4.4804, -4.4606, -4.4484, -4.4241, -4.412,
-1.4018, -1.3949, -1.3582)), class = "data.frame", row.names = c(NA,
-8L))
MGI_Pallador_irri <- read.csv("MGI_Pallador_irri_0.65_late_sowing.csv")
dput(MGI_Pallador_irri)
structure(list(wlieu = c(" 73_063_v3test",
" 73_063_v3test", " 73_063_v3test",
" 73_063_v3test", " 73_063_v3test",
" 73_063_v3test", " 73_063_v3test",
" 73_063_v3test", " 73_063_v3test",
" 73_063_v3test", " 73_063_v3test",
" 73_063_v3test", " 73_063_v3test",
" 73_063_v3test", " 73_063_v3test",
" 73_063_v3test", " 73_063_v3test",
" 73_063_v3test", " 73_063_v3test",
" 73_063_v3test", " 73_073_v3test",
" 73_073_v3test", " 73_073_v3test",
" 73_073_v3test", " 73_073_v3test",
" 73_073_v3test", " 73_073_v3test",
" 73_073_v3test", " 73_073_v3test",
" 73_073_v3test", " 73_073_v3test",
" 73_073_v3test"), mafruit = c(1.367,
1.626, 1.661, 1.558, 1.145, 1.772, 1.86, 1.776, 2.748, 2.832,
2.705, 2.994, 0.689, 0.846, 1.158, 0.656, 0.699, 0.846, 1.158,
0.656, 2.532, 2.677, 2.569, 2.994, 2.554, 2.775, 2.591, 2.94,
2.468, 2.835, 2.665, 2.818)), class = "data.frame", row.names = c(NA,
-32L))
>
MGI_Pallador_no_irri <- read.csv("MGI_Pallador_no_irri_0.65_late_sowing.csv")
dput(MGI_Pallador_no_irri)
structure(list(wlieu = c(" 73_063_v3test",
" 73_063_v3test", " 73_063_v3test",
" 73_063_v3test", " 73_063_v3test",
" 73_063_v3test", " 73_063_v3test",
" 73_063_v3test", " 73_063_v3test",
" 73_063_v3test", " 73_063_v3test",
" 73_063_v3test", " 73_063_v3test",
" 73_063_v3test", " 73_063_v3test",
" 73_063_v3test", " 73_063_v3test",
" 73_063_v3test", " 73_063_v3test",
" 73_063_v3test", " 73_073_v3test",
" 73_073_v3test", " 73_073_v3test",
" 73_073_v3test", " 73_073_v3test",
" 73_073_v3test", " 73_073_v3test",
" 73_073_v3test", " 73_073_v3test",
" 73_073_v3test", " 73_073_v3test",
" 73_073_v3test"), mafruit = c(1.098,
1.235, 1.582, 0.981, 1.145, 1.772, 1.86, 1.776, 1.664, 1.981,
2.193, 2.674, 0.689, 0.846, 1.158, 0.656, 0.699, 0.846, 1.158,
0.656, 2.068, 2.144, 1.97, 2.994, 1.599, 1.693, 1.216, 2.3, 2.108,
2.734, 2.294, 2.726)), class = "data.frame", row.names = c(NA,
-32L))
Solution
First lets break this down into some functions, the first one takes an input file path and returns a raster:
get_simulation <- function(input_file) {
Results <- read.csv(file = input_file, sep = ",", dec = ".", header = TRUE)
Results$CLocation <- gsub(" ", "", Results$wlieu) # 去除列 wlieu 的空格字符
medianYield <- aggregate(mafruit ~ CLocation, data = Results, FUN = median)
[...plus more lines from your source up to...]
raster_Results_wgs84 <- projectRaster(raster_Results, crs = "+proj=longlat +datum=WGS84")
return(raster_Results_wgs84)
}
Then get France with:
get_france = function( ){
url <- "https://raw.githubusercontent.com/gregoiredavid/france-geojson/master/regions.geojson"
france <- st_read(url)
return(france) # this is in wgs84
}
france = get_france()
Now read the files and get the rasters:
r_irri = get_simulation(files[1])
r_nonirri = get_simulation(files[2])
As supplied, one raster is bigger than the other everywhere, so for testing lets tweak one cell:
r_nonirri[5,1] = 2
Now convert the difference to a data frame, remove NAs, and create a factor column which is what we're going to use in ggplot:
diffs = as.data.frame(r_irri > r_nonirri, na.rm=TRUE, xy=TRUE)
diffs$Status = factor(c("Irrigated","Non Irrigated")[diffs$layer+1])
Then this ggplot get the fundamentals:
ggplot() + geom_sf(data=france) +
geom_tile(data=diffs, aes(x=x,y=y,fill=Status))
Garnish that with whatever other ggplot herbs and spices you need.
Not sure what trouble you are having with tmap, but a basic example is:
tm_shape(france) + tm_polygons() +
tm_shape(r1>r2) + tm_raster(labels=c("Non-Irri","Irri"))
which can be styled and coloured to taste. But maybe further tmap problems warrant a new question (and maybe on http://gis.stackexchange.com)
- Restrict R functions for an exam
- How to pass unquoted argument to filter() within user defined function
- Adding a row at the end of a group in R
- How to find the percentile for each observation in a data frame in R?
- str_view() in pdf slides is no longer working
- Reading C style date and time from a binary file using R
- How to create some of the exotic types?
- Making list of strings while summarizing with dplyr
- Randomly fill an zero matrix with some ones with row and column constraints
- Do not understand error message about namespace R with spatialEco package
- Imputing date in time series dataframe
- How to rbind all the data.frames in your working environment?
- Grepl with logical operator AND across multiple alphanumeric columns
- Improving the R traceback when rendering R Markdown non-interactively
- How do I install an R package from the source tarball on windows?
- Create an R function which dinamycally and sequentially applies other predefined functions to an initial input variable
- Percentage labels per x-axis group in stacked bar plots
- Extract rows for the first week of data that meets criteria in R
- R fastest bivariate regression slope coefficient
- Different x-ticks per facet in ggplot2
- Filling in a column in a dataframe using another dataframe that partially matches
- R: Hashing rows of data.table
- How to save changes in a subset of data in a Shiny app with editable DataTable?
- conditionally debug inside a lapply loop in R
- Syntax error in R code executed with Slurm and Singularity
- Header font and merge vertically in flextable in R
- Basic return behavior of reactive
- Optimization of tail recursion in R
- Adding space between bars in ggplot2
- How to Nest Data for Pairwise Group Comparisons in dplyr?