使用栅格计算农药的平均施用量和总施用量，但未将数字相加

如何解决使用栅格计算农药的平均施用量和总施用量，但未将数字相加

我使用一个栅格文件，该文件具有2015年草甘膦对大豆的全球施用量（以千克/公顷为单位）。我想计算每个国家的平均施用量，并得出以千克为单位的总施用量。

我尝试提取数据，但是当我检查数据时，它们并没有加起来，因此我需要一些错误信息。

这是我的代码：

## 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,fun = sum,na.rm = TRUE)

# Get the total area 
total_area_extract = raster::extract(area(rast_trimmed),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）。

但是问题还在于整个应用程序似乎至少具有一个数量级。根据{{​​3}}，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,fun = sum,na.rm = TRUE)

我认为您在这里犯了一个错误。您需要使用非仅NA的单元格。

rarea <- mask(a,r)
tot_area <- raster::extract(rarea,na.rm = TRUE)
## not 
## tot_area <- raster::extract(area(r),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