---

■ Python GIS: Interpolating and Plotting Corn Grain Yield Data

---

In my previous post, I explained how to create a GIS map using Python. Today, I’ll introduce how to create the same GIS map using R.

First, let’s install all the required packages.

if(!require(ggplot2)) install.packages("ggplot2")
library(ggplot2)
if(!require(akima)) install.packages("akima")
library(akima)
if(!require(RCurl)) install.packages("RCurl")
library(RCurl)

and I’ll upload a dataset for practice.

if(!require(readr)) install.packages("readr")
library(readr)
github="https://raw.githubusercontent.com/agronomy4future/raw_data_practice/main/grain_yield_map.csv"
df= data.frame(read_csv(url(github), show_col_types=FALSE))

head (df,5)
  Latitude Longitude     GY
1 12.15725 -106.1403 16.248
2 12.15724 -106.1399 26.703
3 12.15724 -106.1395 16.569
4 12.15723 -106.1391 21.451
5 12.15722 -106.1387 19.107
.
.
.

Next, I’ll extract columns for latitude, longitude, and y (output)

latitude= df$Latitude
longitude= df$Longitude
output= df$GY

and I’ll interpolate data

grid_data=interp(longitude, latitude, output, nx=70, ny=70)
grid_df=expand.grid(x=grid_data$x, y=grid_data$y)
grid_df$z=as.vector(grid_data$z)
grid_df=na.omit(grid_df)

Finally, I’ll create a GIS map using ggplot().

ggplot() +
  geom_raster(data=grid_df, aes(x=x, y=y, fill=z)) + 
  scale_fill_gradientn(colors=c("darkred", "yellow", "darkgreen", "green"), 
                       name="Corn Grain Yield (Mg/ha)", 
                       limits=c(0, 50),
                       breaks=seq(0, 50, 10)) +
  geom_point(data=df, aes(x=Longitude, y=Latitude), color="black", size=2) +
  theme_void() +
windows(width=10, height=5)

---

Full code

If you copy and paste the full code below into your R console, you can obtain the same GIS map.

if(!require(ggplot2)) install.packages("ggplot2")
library(ggplot2)
if(!require(akima)) install.packages("akima")
library(akima)
if(!require(RCurl)) install.packages("RCurl")
library(RCurl)

if(!require(readr)) install.packages("readr")
library(readr)
github="https://raw.githubusercontent.com/agronomy4future/raw_data_practice/main/grain_yield_map.csv"
df= data.frame(read_csv(url(github), show_col_types=FALSE))

# Extract columns for latitude, longitude, and y (output)
latitude= df$Latitude
longitude= df$Longitude
output= df$GY

# Interpolate data
grid_data=interp(longitude, latitude, output, nx=70, ny=70)
grid_df=expand.grid(x=grid_data$x, y=grid_data$y)
grid_df$z=as.vector(grid_data$z)
grid_df=na.omit(grid_df)

# Plot using ggplot2
ggplot() +
  geom_raster(data=grid_df, aes(x=x, y=y, fill=z)) + 
  scale_fill_gradientn(colors=c("darkred", "yellow", "darkgreen", "green"), 
                       name="Corn Grain Yield (Mg/ha)", 
                       limits=c(0, 50),
                       breaks=seq(0, 50, 10)) +
  geom_point(data=df, aes(x=Longitude, y=Latitude), color="black", size=2) +
  theme_void() +
windows(width=10, height=5)

https://github.com/agronomy4future/r_code/blob/main/R_GIS_Interpolating_and_Plotting_Corn_Grain_Yield_Data.ipynb

---

---