之前给大家分享了「1940年～2024年各省市区县10米高度的风速、边界层高度及通风系数月度与年度面板数据」：
1940年～2024年各省市区县10米高度的风速、边界层高度及通风系数月度与年度面板数据: https://rstata.duanshu.com/#/brief/course/5e8f37bf5fc44d329a6efaa466e939ac
今天我们再分享下该数据的下载和处理方法。
通风系数又称空气流通系数，是环境规制研究常用的工具变量。例如：
沈坤荣, 金刚, 方娴. 环境规制引起了污染就近转移吗? 经济研究 2017;52; 44-59.
下面是两篇文章中对通风数据的描述：
通风系数在大气污染的标准方框模型中被认为是空气污染扩散速度的决定因素（Jacobsen，2002）。这个系数被定义为风速和混合高度的乘积，前者决定了污染的水平扩散，后者决定了污染物在大气中扩散的高度。对于两个排放水平相同的地方，通风系数高的地方受空气污染的影响较小。另外由于通风系数是由大规模的天气系统决定的，可以被认为是当地经济活动的外生因素。我们将城市间空气污染差异的这一外生源作为 TCZ 状态的工具。
Broner等（2012）使用的 ERA-Interim 数据提供了75°×75°单元（约83平方公里）的全球网格 10 米高度的风速和混合高度。每个网格单元的通风系数是由平均风速和边界层高度相乘而得到的。然后我们对 1991 年至 1996 年（政策实施前两年）按单元平均这一指标。ERA-Interim 数据库中的位置可以通过经纬度与我们的中国城市相匹配（从 world-gazetteer.com 获得）。我们将每个城市的通风系数定义为其在 ERA-Interim 网格中最接近的四个单元的平均通风系数。
该数据的处理方法也非常简单，首先从 https://cds.climate.copernicus.eu/datasets/reanalysis-era5-single-levels-monthly-means 下载 ERA5 月度数据，只需要其中的 10 米高度的风速（10m Windspeed）、边界层高度（Boundary layer height）两个指标。下载得到的数据是 GRIB 文件（也可以下载 netCDF 格式的 nc 文件），可以作为栅格数据使用 R 语言处理，分别对两个指标进行分省市区县平均再相乘即可得到通风系数指标。
年度通风系数可以通过平均月度通风系数得到。
下载 该数据的原始数据是这里下载的：ERA5 monthly averaged data on single levels from 1940 to present：https://cds.climate.copernicus.eu/datasets/reanalysis-era5-single-levels-monthly-means
首先注册登录，点击 Download 表单：
然后 Product type 里面选择：Monthly averaged reanalysis；Variable 里面选择 Wind 里面的 10m wind speed 和 Other 里面的 Boundary layer height；Year 全选或者只选择自己需要的；Month 全选或只选择自己需要的；Time 选 00:00；Geographical area 可以选全部，也可以选择中国范围的，如果选择中国范围的话，中国的范围是这样的：
library( tidyverse) library( sf) read_sf( "2021行政区划/省.shp" ) %>% st_transform( 4326 ) -> prov prov %>% st_union( ) %>% nngeo:: st_remove_holes( ) %>% st_sf( ) -> cnpoly st_bbox( cnpoly)
#> xmin ymin xmax ymax #> 73.50235 3.83703 135.09567 53.56362
Data format 选择两个都行，都可以处理，这里我选择的是 GRIB；Download format 选择 zip，然后就可以下载了。
下载解压之后就会得到一个 data.grib 文件。
处理 grid 文件可以使用 terra 包的 rast() 函数读取为一个栅格数据：
library( tidyverse) library( sf) library( terra) rast( "data.grib" ) -> rst rst
#> class : SpatRaster #> size : 199, 247, 2054 (nrow, ncol, nlyr) #> resolution : 0.2500041, 0.2500051 (x, y) #> extent : 73.377, 135.128, 3.711997, 53.463 (xmin, xmax, ymin, ymax) #> coord. ref. : lon/lat Coordinate System imported from GRIB file #> source : data.grib #> names : SFC (~[m/s], SFC (~t [m], SFC (~[m/s], SFC (~t [m], SFC (~[m/s], SFC (~t [m], ... #> unit : m/s, m, m/s, m, m/s, m, ... #> time : 1940-01-01 to 2025-07-01 UTC (1027 steps)
可以看到该栅格数据又 time 属性，可以使用 time() 函数提取：
time( rst) %>% as.character ( ) -> timevct length ( timevct)
每个月份都有两个，分别是选择的两个变量：
names ( rst) %>% str_extract( "(Boundary layer height)|(10m\\. Windspeed)" ) %>% str_remove_all( "\\." ) %>% map_chr( ~ str_replace_all( .x, " " , "_" ) ) %>% map_chr( ~ paste0( "sfc_" , .x) ) -> attrvct length ( attrvct)
可以把时间 + 变量名称作为栅格数据的 layer 名称：
rst %>% `names<-`( paste0( attrvct, "__" , timevct) ) -> rst rst
#> class : SpatRaster #> size : 199, 247, 2054 (nrow, ncol, nlyr) #> resolution : 0.2500041, 0.2500051 (x, y) #> extent : 73.377, 135.128, 3.711997, 53.463 (xmin, xmax, ymin, ymax) #> coord. ref. : lon/lat Coordinate System imported from GRIB file #> source : data.grib #> names : sfc_1~01-01, sfc_B~01-01, sfc_1~02-01, sfc_B~02-01, sfc_1~03-01, sfc_B~03-01, ... #> unit : m/s, m, m/s, m, m/s, m, ... #> time : 1940-01-01 to 2025-07-01 UTC (1027 steps)
然后就可以对栅格数据进行省市区县的汇总了：
读取省市区县矢量数据：
合并省份的数据可以得到中国的范围了，进而提取中国区域的栅格数据保存：
prov %>% st_union( ) %>% nngeo:: st_remove_holes( ) -> cn rst %>% terra:: crop( vect( cn) ) %>% terra:: mask( vect( cn) ) -> cnrst cnrst %>% writeRaster( "cnrst.tif" , overwrite = TRUE )
使用 terra::extract() 分区域求均值：
处理提取结果。在 ERA-Interim 数据中，Boundary_layer_height（边界层高度）指的是大气边界层（Atmospheric Boundary Layer, ABL）的顶部高度。大气边界层是靠近地表的一层大气，受地表摩擦、加热和冷却等过程直接影响。把 10m 高风速和边界层高度相乘就得到了通风系数（空气流通系数）。
provres %>% left_join( prov %>% st_drop_geometry( ) %>% mutate( ID = row_number( ) ) ) %>% select( - ID) %>% select( 省, 省代码, everything( ) ) %>% gather( 3 : ncol( .) , key = "key" , value = "value" ) %>% separate( key, into = c ( "class" , "日期" ) , sep = "__" ) %>% mutate( 日期 = ymd( 日期) , 年份 = year( 日期) ) %>% spread( class , value) %>% mutate( ventilation_coefficient = sfc_10m_Windspeed * sfc_Boundary_layer_height) -> provdf provdf
#> # A tibble: 34,918 × 7 #> 省 省代码 日期 年份 sfc_10m_Windspeed sfc_Boundary_layer_height #> <chr> <dbl> <date> <dbl> <dbl> <dbl> #> 1 上海市 310000 1940-01-01 1940 4.02 584. #> 2 上海市 310000 1940-02-01 1940 3.91 497. #> 3 上海市 310000 1940-03-01 1940 3.44 481. #> 4 上海市 310000 1940-04-01 1940 3.90 489. #> 5 上海市 310000 1940-05-01 1940 3.86 462. #> 6 上海市 310000 1940-06-01 1940 3.97 441. #> 7 上海市 310000 1940-07-01 1940 5.08 587. #> 8 上海市 310000 1940-08-01 1940 3.65 539. #> 9 上海市 310000 1940-09-01 1940 4.32 656. #> 10 上海市 310000 1940-10-01 1940 3.21 534. #> # ℹ 34,908 more rows #> # ℹ 1 more variable: ventilation_coefficient <dbl>
cityres %>% left_join( city %>% st_drop_geometry( ) %>% mutate( ID = row_number( ) ) ) %>% select( - ID) %>% select( 省, 省代码, 市, 市代码, everything( ) ) %>% gather( 5 : ncol( .) , key = "key" , value = "value" ) %>% separate( key, into = c ( "class" , "日期" ) , sep = "__" ) %>% mutate( 日期 = ymd( 日期) , 年份 = year( 日期) ) %>% spread( class , value) %>% mutate( ventilation_coefficient = sfc_10m_Windspeed * sfc_Boundary_layer_height) -> citydf citydf
#> # A tibble: 381,017 × 9 #> 省 省代码 市 市代码 日期 年份 sfc_10m_Windspeed #> <chr> <dbl> <chr> <dbl> <date> <dbl> <dbl> #> 1 安徽省 340000 安庆市 340800 1940-01-01 1940 2.95 #> 2 安徽省 340000 蚌埠市 340300 1940-01-01 1940 3.11 #> 3 安徽省 340000 亳州市 341600 1940-01-01 1940 3.00 #> 4 安徽省 340000 池州市 341700 1940-01-01 1940 2.66 #> 5 安徽省 340000 滁州市 341100 1940-01-01 1940 3.15 #> 6 安徽省 340000 阜阳市 341200 1940-01-01 1940 3.12 #> 7 安徽省 340000 合肥市 340100 1940-01-01 1940 3.16 #> 8 安徽省 340000 淮北市 340600 1940-01-01 1940 2.95 #> 9 安徽省 340000 淮南市 340400 1940-01-01 1940 3.08 #> 10 安徽省 340000 黄山市 341000 1940-01-01 1940 2.13 #> # ℹ 381,007 more rows #> # ℹ 2 more variables: sfc_Boundary_layer_height <dbl>, #> # ventilation_coefficient <dbl>
countyres %>% left_join( county %>% st_drop_geometry( ) %>% mutate( ID = row_number( ) ) ) %>% select( - ID) %>% select( 省, 省代码, 市, 市代码, 县, 县代码, everything( ) ) %>% gather( 7 : ncol( .) , key = "key" , value = "value" ) %>% separate( key, into = c ( "class" , "日期" ) , sep = "__" ) %>% mutate( 日期 = ymd( 日期) , 年份 = year( 日期) ) %>% spread( class , value) %>% mutate( ventilation_coefficient = sfc_10m_Windspeed * sfc_Boundary_layer_height) -> countydf countydf
#> # A tibble: 2,954,679 × 11 #> 省 省代码 市 市代码 县 县代码 日期 年份 sfc_10m_Windspeed #> <chr> <dbl> <chr> <dbl> <chr> <dbl> <date> <dbl> <dbl> #> 1 安徽省 340000 安庆市 340800 大观区 340803 1940-01-01 1940 3.49 #> 2 安徽省 340000 安庆市 340800 怀宁县 340822 1940-01-01 1940 3.38 #> 3 安徽省 340000 安庆市 340800 潜山市 340882 1940-01-01 1940 2.68 #> 4 安徽省 340000 安庆市 340800 宿松县 340826 1940-01-01 1940 3.40 #> 5 安徽省 340000 安庆市 340800 太湖县 340825 1940-01-01 1940 2.58 #> 6 安徽省 340000 安庆市 340800 桐城市 340881 1940-01-01 1940 2.95 #> 7 安徽省 340000 安庆市 340800 望江县 340827 1940-01-01 1940 3.57 #> 8 安徽省 340000 安庆市 340800 宜秀区 340811 1940-01-01 1940 3.49 #> 9 安徽省 340000 安庆市 340800 迎江区 340802 1940-01-01 1940 3.40 #> 10 安徽省 340000 安庆市 340800 岳西县 340828 1940-01-01 1940 2.06 #> # ℹ 2,954,669 more rows #> # ℹ 2 more variables: sfc_Boundary_layer_height <dbl>, #> # ventilation_coefficient <dbl>
最后保存数据：
provdf %>% haven:: write_dta( "1940年1月～2025年7月各省份10米高度的风速、边界层高度及通风系数面板数据.dta" , label = "数据计算：微信公众号 RStata" ) citydf %>% haven:: write_dta( "1940年1月～2025年7月各城市10米高度的风速、边界层高度及通风系数面板数据.dta" , label = "数据计算：微信公众号 RStata" ) countydf %>% haven:: write_dta( "1940年1月～2025年7月各区县10米高度的风速、边界层高度及通风系数面板数据.dta" , label = "数据计算：微信公众号 RStata" )
年度数据可以平均月度数据得到：
这样我们就处理得到了年度和月度的数据。
点击这里跳转到 RStata 短书平台获取附件：如何使用 R 语言计算各省市区县空气流通系数
评论