在论文「超时加班与劳动收入份额：基于卫星夜间灯光的经验证据」中，作者提出了这么一种判断企业是否加班的方法：
简单来说就是有三种标准：
- 时间维度标准：该公司处的夜间灯光亮度大于该年节假日夜间灯光亮度的中位数；
- 空间维度标准：该公司处的夜间灯光亮度大于周边8个网格的中位数；
- 时空双维度：同时满足上面的时间和空间标准的日子才被视为加班日。
今天的课程中我们将会讲解如何使用 R 语言提取上市公司所在位置及周边 8 个格点的夜间灯光亮度、判断是否加班及统计年加班天数。
在之前的课程中我们讲解过如何使用 R 语言爬取和处理 VNP46A2 的日度夜间灯光亮度数据：
名师讲堂｜使用 R 语言下载和处理日度夜间灯光栅格数据：https://rstata.duanshu.com/#/brief/course/270178a4d6bb4ee2ba8b66ae71de3b8f
基于该课程讲解的方法处理得到的栅格数据可以从这里下载到：
2012～2024 年 VNP46A2 日度夜间灯光亮度栅格数据：https://rstata.duanshu.com/#/brief/course/ae5198d142b944108bcdecbd56e1fe36
注意：该夜间灯光亮度的单位是 0.1 nWatts/(cm^2 sr) 如果想转换成 nWatts/(cm^2 sr)，需要把数值乘以 0.1 。
今天我们将会在该课程的基础上讲解。
除了日度夜间灯光栅格数据外，还会需要下面两个数据：
- 2000～2023 年上市公司注册地址与办公地址（含经纬度及其所处的省市区县）：https://rstata.duanshu.com/#/brief/course/90fa75897e564a45beb32f7840be0787
- 2007~2024 年法定节假日及星期数据：https://rstata.duanshu.com/#/brief/course/63147c1f06924ddb89b152d1df84b78c
附件中我也给大家提供了示例数据：
- 2024年日度夜间灯光数据（日度夜间灯光栅格数据）
- 2024年法定节假日及星期数据.dta
- 2023年上市公司办公地址经纬度.dta
由于截止讲义材料编写的时候还没有 2024 年的上市公司办公地址数据，所以我们只能假设 2023～2024 年上市公司没有搬家，也就是使用 23 年的数据作为 24 年的数据替代。
提取上市公司所在位置和周边 8 个相邻格点的夜光亮度
首先加载所需 R 包：
library(tidyverse) |
索引所有的日度夜间灯光栅格数据：
fs::dir_ls("2024年日度夜间灯光数据", regexp = "tif$") -> fls |
不过 fs::dir_ls() 函数在 Windows 电脑上使用可能会遇到中文乱码问题，可以通过下面两种方法解决：
- 把文件夹的名字还有文件的名字都换成不带中文的；
- 改用 list.files();
paste0("2024年日度夜间灯光数据/", list.files("2024年日度夜间灯光数据")) -> fls |
这样效果也是一样的，不过麻烦一些。
读取一个栅格数据：
rast(fls[length(fls)]) -> rst1 |
绘图展示：
read_sf("九段线.geojson") -> jdx |
读取上市公司办公地址经纬度数据：
haven::read_dta("2023年上市公司办公地址经纬度.dta") %>% |
由于 terra 包里面并没有提供直接计算周边 8 个格点亮度中位数的方法，所以我们还是得手动编程解决。
一种办法是下面这样：
library(terra) |
# 获取某个特定点的值（例如第 5 行第 5 列的点） |
focal()函数是用来进行邻域计算的，这个方法看起来很方便，不过实际上我也这么做了，结果程序运行了半个月才运行完。
最初我以为
focal(r, w=3, median, na.rm)就可以了，后来发现这是计算是计算九个格点共同中位数的。
由于上述的方法极度低效，所以我们还是得更换方法。最后我发现还是使用矩阵进行提取最快。
这个图和分享的数据里面的不太一样是因为那个图绘制的时候给平安银行的经纬度设置成了早年的。
上图展示了「平安银行」的经纬度及周边的夜光栅格。假设我们的栅格数据就这么大，我们可以数一下，这个栅格的列数为：
- ncol = 6;
如果从左到右对每个格子进行编号，那么可以数一数，平安银行所在格点的编号是：
- 平安银行所在格点的编号是: 16;
那么它周边 8 个格点的编号分别为：
- grid1: 16 - 6 - 1 = 9
- grid2: 16 - 6 = 10
- grid3: 16 - 6 + 1 = 11
- grid4: 16 - 1 = 15
- grid5: 16 = 16
- grid6: 16 + 1 = 17
- grid7: 16 + 6 - 1 = 21
- grid8: 16 + 6 = 22
- grid9: 16 + 6 + 1 = 23
因此我们只需要知道栅格数据的列数以及某个格点的索引就可以推出其周边 8 个格点的索引。
那么在整个栅格数据里面，平安银行所在格点为：
xy = cbind(lon = df$办公地址_经度[1], lat = df$办公地址_纬度[1]) |
# 栅格数据转换成一维的矩阵 |
# 该处的值为 |
其周边 8 个格点的索引分别为：
# 栅格数据的列是 |
所以所有经纬度及周边网格的索引是：
df %>% |
提取所有上市公司该日 9 个格点的夜光亮度：
df2 %>% |
然后我们就可以循环所有的日子了，为了更高效，这里我用的是多线程循环：
dir.create("res") |
合并提取结果：
fs::dir_ls("res") %>% |
计算中位数与比较是否加班
计算周边 8 个网格的中位数：
df %>% |
上市公司所在位置的：
df %>% |
读取节假日数据：
haven::read_dta("2024年法定节假日及星期数据.dta") %>% |
合并上面的三个数据：
df1 %>% |
空间维度上的比较
dfall %>% |
时间维度上的比较
# 计算年度法定节假日灯光亮度中位数 |
同时满足两个维度
dfall3 %>% |
如果想统计加班天数，分年汇总即可：
dfall3 %>% |
如果想计算加班天数的比例，再计算一个各年的天数即可：
dfall3 %>% |
这样我们就计算得到了这个指标～
最后再补充下上面图表的绘制代码：
library(tidyverse) |
如何参加课程？
购买 RStata 名师讲堂会员即可参加该课程啦（之前的和未来的都可以参加）！
价格：2800/年 或者 4800/长期
购买会员可以从这里下单：https://rstata.duanshu.com/#/card/list/
名师讲堂会员权益：
- 参加每个月 3～4 次的名师讲堂课程；
- 参加平台上的其他 R 语言和 Stata 的课程；
- 以会员折扣价购买我们分享的数据资料（10 元/份）；
- 课程内外的提问解答服务（课程外的尽量帮忙解决）。
* 如果发票可添加小编微信 r_stata （RStata 李老师）开具。如需数据资料，购买后可添加小编微信免费领取数据折扣卡。
更多关于 RStata 会员的更多信息可添加微信号 r_stata 咨询：
课程主页（点击文末的阅读原文即可跳转）：
点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 R 语言提取上市公司所在位置及周边 8 个格点的夜间灯光亮度、判断是否加班及统计年加班天数
评论