用栅格计算农药的平均施用量和总施用量，但数字不相加
Using raster to calculate the mean application and total application of pesticides, but numbers not adding up
我使用了一个栅格文件，其中包含 2015 年草甘膦对大豆的全球施用率（在 kg/hectare 中）。我想计算每个国家的平均施用率，并获得总施用量公斤
我试图提取数据，但当我检查它时，它并没有加起来，所以我需要一些帮助来弄清楚我哪里出错了。
这是我的代码：
## Load libraries and read data ----
library(raster)
library(tidyverse)
library(maptools) # To get the borders of countries
# Border of countries
data(wrld_simpl)
# Make a data frame to put my numbers on the corresponding countries
world_data_empty = data.frame(country = wrld_simpl$NAME, iso3c = wrld_simpl$ISO3)
# This is a file with the application rate of glyphosate to soybeans in 2015 in kg/ha
raster_path = "https://raw.github.com/hansronald/Pesticide-data/master/APR_Soybean_Glyphosate_2015_L.tif"
sqm_to_ha_conversion = 0.0001
rast <- raster(raster_path)
## Extract the raster data ----
# Remove all the negative values of the raster
# This is because negative values correspond to things like water and I dont want them to count when adding, replace with NAs
# Then trim to remove NAs
rast_positive = rast
rast_positive <- clamp(rast, lower=0, useValues=FALSE)
rast_trimmed = trim(rast_positive)
# Multiply the kg/ha with the area of each cell to get the total kgs
rast_total_pesticide_application = rast_trimmed * area(rast_trimmed)
# Get the mean applcation rate (kg/ha) for each country
mean_application_rate_extract = raster::extract(rast_trimmed, wrld_simpl, fun = mean, na.rm = TRUE) # sp = T for keeping original dataframe
# Get the total applcation rate (kg/ha) for each country
total_application_extract = raster::extract(rast_total_pesticide_application, wrld_simpl, fun = sum, na.rm = TRUE)
# Get the total area
total_area_extract = raster::extract(area(rast_trimmed), wrld_simpl, fun = sum, na.rm = TRUE)
## Create the data frame ----
# Put the mean pesticide use per country in a dataframe and name the column after the pesticide
# NaNs created because all raster cells are NA in country polygon, replace with 0
mean_application_rate_extract[is.nan(mean_application_rate_extract)] = 0
mean_application_rate = data.frame(mean_application_rate = mean_application_rate_extract)
total_application = data.frame(total_application = total_application_extract)
total_area = data.frame(total_area = total_area_extract)
world_data = data.frame(world_data_empty, mean_application_rate, total_application, total_area) %>%
as_tibble()
# Check calculations by selecting a few countries and multiplying apr*area
world_data %>%
filter(iso3c %in% c("CHN", "BRA", "USA")) %>%
mutate(total_application_calc = mean_application_rate * total_area)
这是输出
country iso3c mean_application_rate total_application total_area total_application_calc
<fct> <fct> <dbl> <dbl> <dbl> <dbl>
1 Brazil BRA 2.00 4253187. 84469. 168653.
2 China CHN 2.09 9153007. 93254. 194736.
3 United States USA 1.93 5070446. 93889. 181164.
所以这里有一些问题。首先是 total_application_calc 应该等于总应用（因为它是应用率 (kg/ha) 乘以总面积 (ha)。
但问题还在于，整个应用程序似乎至少有一个数量级。根据 this data，2014 年草甘膦在大豆上的总应用量为 122,473,987 磅，相当于 55,553,266 公斤，而我从该数据集中得到的是 5,070,446 公斤。如果它稍微偏离也没关系，因为它们是不同的来源和不同的假设，但不是那么多。
任何人都可以帮助我做错了什么吗？
我稍微简化了您的代码，我认为这些数字现在更有意义了。我无法回答外部有效性问题。
library(raster)
library(maptools)
data(wrld_simpl)
r <- raster("https://raw.github.com/hansronald/Pesticide-data/master/APR_Soybean_Glyphosate_2015_L.tif")
r <- clamp(r, lower=0, useValues=FALSE)
# area in ha
a <- area(r) * 100
mean_app <- raster::extract(r, wrld_simpl, fun = mean, na.rm = TRUE)
rtot <- r * a
tot_app <- raster::extract(rtot, wrld_simpl, fun = sum, na.rm = TRUE)
我认为你在这里犯了一个错误。您需要使用不仅仅是 NA 的单元格。
rarea <- mask(a, r)
tot_area <- raster::extract(rarea, wrld_simpl, fun = sum, na.rm = TRUE)
## not
## tot_area <- raster::extract(area(r), wrld_simpl, fun = sum, na.rm = TRUE)
检查结果
w <- data.frame(country = wrld_simpl$NAME, iso3c = wrld_simpl$ISO3,
mean_app=mean_app, tot_app=tot_app, tot_area=tot_area)
w$tot_calc <- w$tot_area * w$mean_app
w[is.na(w)] <- 0
w[w$iso3c %in% c("CHN", "BRA", "USA"), ]
# country iso3c mean_app tot_app tot_area tot_calc
#21 Brazil BRA 1.996631 425318662 213982181 427243427
#30 China CHN 2.088219 915300667 439036703 916804725
#209 United States USA 1.929559 507044556 263544939 508525402
global_app <- cellStats(rtot, "sum")
global_app
# 2375398749
sum(w$tot_app)
# 2367120358
我使用了一个栅格文件，其中包含 2015 年草甘膦对大豆的全球施用率（在 kg/hectare 中）。我想计算每个国家的平均施用率，并获得总施用量公斤
我试图提取数据，但当我检查它时，它并没有加起来，所以我需要一些帮助来弄清楚我哪里出错了。
这是我的代码：
## Load libraries and read data ----
library(raster)
library(tidyverse)
library(maptools) # To get the borders of countries
# Border of countries
data(wrld_simpl)
# Make a data frame to put my numbers on the corresponding countries
world_data_empty = data.frame(country = wrld_simpl$NAME, iso3c = wrld_simpl$ISO3)
# This is a file with the application rate of glyphosate to soybeans in 2015 in kg/ha
raster_path = "https://raw.github.com/hansronald/Pesticide-data/master/APR_Soybean_Glyphosate_2015_L.tif"
sqm_to_ha_conversion = 0.0001
rast <- raster(raster_path)
## Extract the raster data ----
# Remove all the negative values of the raster
# This is because negative values correspond to things like water and I dont want them to count when adding, replace with NAs
# Then trim to remove NAs
rast_positive = rast
rast_positive <- clamp(rast, lower=0, useValues=FALSE)
rast_trimmed = trim(rast_positive)
# Multiply the kg/ha with the area of each cell to get the total kgs
rast_total_pesticide_application = rast_trimmed * area(rast_trimmed)
# Get the mean applcation rate (kg/ha) for each country
mean_application_rate_extract = raster::extract(rast_trimmed, wrld_simpl, fun = mean, na.rm = TRUE) # sp = T for keeping original dataframe
# Get the total applcation rate (kg/ha) for each country
total_application_extract = raster::extract(rast_total_pesticide_application, wrld_simpl, fun = sum, na.rm = TRUE)
# Get the total area
total_area_extract = raster::extract(area(rast_trimmed), wrld_simpl, fun = sum, na.rm = TRUE)
## Create the data frame ----
# Put the mean pesticide use per country in a dataframe and name the column after the pesticide
# NaNs created because all raster cells are NA in country polygon, replace with 0
mean_application_rate_extract[is.nan(mean_application_rate_extract)] = 0
mean_application_rate = data.frame(mean_application_rate = mean_application_rate_extract)
total_application = data.frame(total_application = total_application_extract)
total_area = data.frame(total_area = total_area_extract)
world_data = data.frame(world_data_empty, mean_application_rate, total_application, total_area) %>%
as_tibble()
# Check calculations by selecting a few countries and multiplying apr*area
world_data %>%
filter(iso3c %in% c("CHN", "BRA", "USA")) %>%
mutate(total_application_calc = mean_application_rate * total_area)
这是输出
country iso3c mean_application_rate total_application total_area total_application_calc
<fct> <fct> <dbl> <dbl> <dbl> <dbl>
1 Brazil BRA 2.00 4253187. 84469. 168653.
2 China CHN 2.09 9153007. 93254. 194736.
3 United States USA 1.93 5070446. 93889. 181164.
所以这里有一些问题。首先是 total_application_calc 应该等于总应用（因为它是应用率 (kg/ha) 乘以总面积 (ha)。
但问题还在于，整个应用程序似乎至少有一个数量级。根据 this data，2014 年草甘膦在大豆上的总应用量为 122,473,987 磅，相当于 55,553,266 公斤，而我从该数据集中得到的是 5,070,446 公斤。如果它稍微偏离也没关系，因为它们是不同的来源和不同的假设，但不是那么多。
任何人都可以帮助我做错了什么吗？
我稍微简化了您的代码，我认为这些数字现在更有意义了。我无法回答外部有效性问题。
library(raster)
library(maptools)
data(wrld_simpl)
r <- raster("https://raw.github.com/hansronald/Pesticide-data/master/APR_Soybean_Glyphosate_2015_L.tif")
r <- clamp(r, lower=0, useValues=FALSE)
# area in ha
a <- area(r) * 100
mean_app <- raster::extract(r, wrld_simpl, fun = mean, na.rm = TRUE)
rtot <- r * a
tot_app <- raster::extract(rtot, wrld_simpl, fun = sum, na.rm = TRUE)
我认为你在这里犯了一个错误。您需要使用不仅仅是 NA 的单元格。
rarea <- mask(a, r)
tot_area <- raster::extract(rarea, wrld_simpl, fun = sum, na.rm = TRUE)
## not
## tot_area <- raster::extract(area(r), wrld_simpl, fun = sum, na.rm = TRUE)
检查结果
w <- data.frame(country = wrld_simpl$NAME, iso3c = wrld_simpl$ISO3,
mean_app=mean_app, tot_app=tot_app, tot_area=tot_area)
w$tot_calc <- w$tot_area * w$mean_app
w[is.na(w)] <- 0
w[w$iso3c %in% c("CHN", "BRA", "USA"), ]
# country iso3c mean_app tot_app tot_area tot_calc
#21 Brazil BRA 1.996631 425318662 213982181 427243427
#30 China CHN 2.088219 915300667 439036703 916804725
#209 United States USA 1.929559 507044556 263544939 508525402
global_app <- cellStats(rtot, "sum")
global_app
# 2375398749
sum(w$tot_app)
# 2367120358