全国分省、市、县NPP面板数据是如何制作的

净初级生产力（Net primary productivity, NPP）是研究陆地生态系统中物质和能量转换的重要指标，NPP的空间分布与区域气候、植被生长以及人类活动等因素息息相关，其变化能反映植被群落的生产能力，是生态系统功能和结构变化的重要表征。

前几天给大家分享了一套NPP数据：【数据分享】2001-2020年中国分省市县净初级生产力（NPP）面板数据工具变量，在这里给大家介绍一下制作方法。

制作思路

1. 获取NPP栅格数据，在这直接使用的是MODIS年度NPP数据产品
2. 对NPP进行分省、分市、分县分区统计

NPP数据获取

MOD17A3HGF V6产品提供了500m分辨率的全球年度NPP数据产品。年度NPP是由8天净光合作用产品（Net Photosynthesis, PSN）MOD17A2H相加合成的。PSN是总初级生产力（Gross Primary Productivity, GPP）和维持呼吸作用（Maintenance Respiration, MR）相减得到的。

library(rgee)
library(sf)

ee_Initialize(drive = T)
# Define a region of interest with sf
ee_roi <- read_sf("D:/R/CMRrgee/SHP/ChinaSizhi.shp") %>%
  sf_as_ee()

npp = ee$ImageCollection("MODIS/006/MOD17A3HGF")$select("Npp")
#数据预览
npp_composite = npp$
  filter(ee$Filter$date('2020-01-01', '2020-12-31'))$
  max()$clip(ee_roi)

nppParams <- list(min = 0, max= 19000, palette = c(
  "#d73027", "#f46d43", "#fdae61",
  "#fee08b", "#d9ef8b", "#a6d96a",
  "#66bd63", "#1a9850"
))

#Map$addLayer(npp_composite, nppParams, "npp")
#Map$centerObject(ee_roi)

cn_npp <- ee_as_raster(
  image = npp_composite,
  region = ee_roi$geometry(),
  scale = 500,
  dsn = "./CNNPP/2020.tif",
  via = "drive",
  maxPixels = 1e+13
)

2020年NPP

上面的代码获取的NPP产品相关参数：

• NPP数据单位：
• 数据尺度：10000

分区统计

由于NPP是分年的时间序列，为了方便获取时间序列面板数据，在这我使用R语言进行分区计算，直接获取2001-2020年的面板数据。

在使用R语言分区计算之前，需要准备分区矢量。

分区矢量准备

由于用到R语言进行分区统计，那么最好是给矢量以英文命名，文件名、字段名都以英文表示。在这里我使用了数读城事提供的2020年行政区划矢量：有需要的请前往这个推文，找菌老板获取【数据分享】2020年度行政区划调整的数据更新与分享（省市县最新）

先把菌老板那的原始数据导出一下，作为实验数据，文件名改为英文。

导出数据

把字段也改为英文，避免在R语言中使用terra包读取时出现报错。

新建字段，复制内容

删除中文字段

NPP分省市县计算

接下来就可以用R语言进行分区统计了。

library(terra)
npplist = list.files(path = "./npp_cn/", pattern = ".tif$")
nppdir = paste0("./npp_cn/", npplist)
npp_cn = rast(nppdir)
names(npp_cn) = seq(2001, 2020, 1)

#分省NPP计算
ShengVect = vect("./SHP/Sheng2020.shp")
ShengRast = rasterize(ShengVect, npp_cn, field="CODE")
NPPSheng = zonal(npp_cn, ShengRast, fun = "mean", na.rm=TRUE)

#分市NPP计算
ShiVect = vect("./SHP/Shi2020.shp")
ShiRast = rasterize(ShiVect, npp_cn, field="CODE")
NPPShi = zonal(npp_cn, ShiRast, fun = "mean", na.rm=TRUE)

#分县NPP计算
XianVect = vect("./SHP/Xian2020.shp")
XianRast = rasterize(XianVect, npp_cn, field="CODE")
NPPXian = zonal(npp_cn, XianRast, fun = "mean", na.rm=TRUE)

#计算结果导出
write.csv(NPPSheng, "./nppdata/NPPSheng.csv")
write.csv(NPPShi, "./nppdata/NPPShi.csv")
write.csv(NPPXian, "./nppdata/NPPXian.csv")

QGIS绘图

分区统计完成的结果，可以使用QGIS挂接一下数据，把分区统计输出的csv数据挂接到SHP上，然后根据NPP数值进行渲染可视化，生成NPP分布图。

QGIS加载csv数据

QGIS数据挂接

绘图的具体操作不再赘述，可以看下面参考文献中QGIS相关推文，或者阅读原文跳转培训班学习。

2020年NPP分县分布

参考文献

1. https://blog.sciencenet.cn/blog-3424049-1264122.html