近期由于需要处理EDGAR网站的nc文件，真是愁坏了我这个没学过地理计算的人。尝试自己根据网上的代码修改，但最终结果一直报错，正巧在统计之都公众号上遇到了这本书的推荐还拿到了签名版，最终根据书中的教程也是顺利解决了各种报错，成功提取了相关数据。如果有需要从EDGAR的nc文件提取区域碳排放的同僚，可以尝试如下R代码：
# nc转tif
# 加载R包
library(rgdal)
library(ncdf4)
library(sp)
library(raster)
library(ggplot2)
year <- seq(2000,2022,1)
year_chart <- as.character(year)
folder_path <- ""
# 遍历文件列表
for (i in year_chart) {
# 读取nc文件
file_path <- paste(folder_path,"EDGAR_2024_GHG_N2O_", i, "_TOTALS_emi.nc", sep="")
# nc转raster
emission = brick(file_path,varname='emissions')
# 计算
emission_sum = calc(emission, fun=sum)
# 导出为tiff
file_out_path <- paste(folder_path,i, sep="")
writeRaster(emission_sum, file.path(file_out_path), format = "GTiff")
}
# tif转excel
# 加载包
library(sp)
library(raster)
library(terra)
library(dplyr)
library(writexl)
library(tidyverse)
library(readxl)
# 设定路径
tif_path <- "X"
shp_path <- "X"
excel_path <- "X"
# 构建循环体
year <- seq(2000,2020,1)
year_char <- as.character(year)
# 以year控制循环
for (i in year_char){
# 加载tif数据
raster_file <- paste(tif_path,i,'.tif',sep="")
raster_data <- rast(raster_file)
# 加载shp数据
shape_file <- paste(shp_path,'China', i ,'County.shp',sep="")
shape_data <- vect(shape_file)
# 投影变换
raster_data_projection <- project(raster_data,crs(shape_data))
# 裁剪
crop_data <- crop(raster_data_projection,shape_data)
# 分区统计
zone_mean <- zonal(crop_data,shape_data,fun='mean')
zone_max <- zonal(crop_data,shape_data,fun='max')
zone_min <- zonal(crop_data,shape_data,fun='min')
zone_sum <- zonal(crop_data,shape_data,fun='sum')
zone_sd <- zonal(crop_data,shape_data,fun='sd')
zone_median <- zonal(crop_data,shape_data,fun='median')
# 聚合及重命名
j <- as.numeric(i)
if(j < 2014){
k <- paste('CODE',i,sep='')
shape_data$YEAR,zone_mean,zone_max,zone_min,zone_sd,zone_sum,zone_median)
column_names <- names(combined_data)
column_names <- c('县级','县级码','year','mean','max','min','sd','sum','median')
names(combined_data) <- column_names
combined_data$year <- j
} else {
combined_data <- cbind(shape_data$县级,shape_data$县级码,
shape_data$year,zone_mean,zone_max,zone_min,zone_sd,zone_sum,zone_median)
column_names <- names(combined_data)
column_names <- c('县级','县级码','year','mean','max','min','sd','sum','median')
names(combined_data) <- column_names
}
# 数据导出
excel_file = paste(excel_path,i,'.xlsx',sep="")
write_xlsx(combined_data,path=excel_file)
}
combined_data <- data.frame()
for (i in year_char){
file_path <- paste(excel_path,i,'.xlsx',sep="")
current_data <- read_excel(file_path,sheet="Sheet1")
combined_data <- bind_rows(combined_data,current_data)
}
write_xlsx(combined_data,path='X.xlsx')
这本书很值得学习，也欢迎大家就上述代码进行交流
