如何根据经纬度判断该地点所处的省份?

如何根据经纬度判断该地点所处的省份?

今天我的琼琼小伙伴给我发了一个数据集,是 2001 年到 2018 年所有上市公司的经纬度数据。她想知道每个公司所处的省份,我就帮她计算了一下。

首先我们先把读取一份中国省份分布的数据,读入之后将坐标系转换为 4326 坐标系(EPSG:4326,+proj=longlat +datum=WGS84 +no_defs)(这个没有必要深究,这种坐标系的特点是经线相互平行)。最后在使用 lwgeom 包的 st_make_valid() 函数矫正 geometry(以防有 invalid geometry):

1
2
3
4
5
6
7
8
library(sf)
library(tidyverse)
library(ggplot2)
library(readxl)
library(lwgeom)
cnmap <- read_sf('100000_full.json') %>%
st_transform(crs = 4326) %>%
st_make_valid()

可以看到 cnmap 对应的投影坐标系是 EPSG 4326:

1
2
3
4
st_crs(cnmap)
#> Coordinate Reference System:
#> EPSG: 4326
#> proj4string: "+proj=longlat +datum=WGS84 +no_defs"

然后我们再把琼琼的这份数据读进来,筛选出 lng 和 lat 变量都不是 NA 的观测值,再使用 st_as_sf() 函数将经纬度坐标变量转变为 geometry 变量(同时这个数据框也就变成了一个 sf 对象,注意这里使用的也是 4326 坐标系),最后使用 st_intersection 将选取所有位于中国境内的观测值:

1
2
3
4
5
df <- read_xlsx('geo.xlsx') %>% 
dplyr::filter(!is.na(lng), !is.na(lat)) %>%
st_as_sf(coords = c("lng", "lat"), crs = 4326) %>%
st_make_valid() %>%
st_intersection(cnmap)

我们画一幅简单的地图来展示这个分布,看了昨天推文的小伙伴估计就注意到了我使用了一个神奇的东西,这是个无法直接显示的字符,这个字符是我在电脑上 Cmd + C 然后 Cmd + V 复制粘贴的,在 FontAwesome 字体中表示地图标记符号:

1
2
3
4
5
6
7
8
9
10
11
12
13
ggplot() + 
geom_sf(data = cnmap, size = 0.2, fill = NA,
color = "white") +
geom_sf_text(data = df,
label = "", family = "FontAwesome",
alpha = 0.4, size = 2,
color = "#E31A1C") +
coord_sf(crs = "+proj=laea +lat_0=23 +lon_0=113 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs") +
theme_modern_rc(base_family = cnfont) +
worldtilegrid::theme_enhance_wtg() +
labs(x = "", y = "",
title = "中国上市公司的分布",
subtitle = "TidyFriday Project")

大家如果无法理解这个图层,可以使用下面的绘图:

1
2
3
4
5
6
7
8
9
10
11
12
ggplot() + 
geom_sf(data = cnmap, size = 0.2, fill = NA,
color = "white") +
geom_sf(data = df, shape = 2,
alpha = 0.4, size = 2,
color = "#E31A1C") +
coord_sf(crs = "+proj=laea +lat_0=23 +lon_0=113 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs") +
theme_modern_rc(base_family = cnfont) +
worldtilegrid::theme_enhance_wtg() +
labs(x = "", y = "",
title = "中国上市公司的分布",
subtitle = "TidyFriday Project")

出来的图是这样的:

下面我们先编写一个简单的函数:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
fun <- function(name){
st_intersection(df, cnmap[cnmap$name == name,]) %>%
st_drop_geometry() %>%
select(code, year)
}

fun("北京市")
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
#> # A tibble: 3,526 x 2
#> code year
#> * <dbl> <dbl>
#> 1 8 2015
#> 2 8 2016
#> 3 8 2017
#> 4 10 2010
#> 5 10 2011
#> 6 10 2013
#> 7 10 2012
#> 8 10 2015
#> 9 10 2014
#> 10 18 2015
#> # … with 3,516 more rows

这个函数的功能是可以输入省级行政单位的名称,然后返回 df 中所有位于这个行政单位内的观测值。st_drop_geometry() 用于删除 sf 对象中的 geometry 变量,同时把这个 sf 对象退化成普通的数据框。

下面我们计算每个观测值的省份:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
library(purrr)
cnmap %>%
select(name, geometry) %>%
mutate(
data = map(name, fun)
) %>%
st_drop_geometry() -> nestdf
nestdf
#> # A tibble: 35 x 2
#> name data
#> * <chr> <list>
#> 1 北京市 <tibble [3,526 × 2]>
#> 2 天津市 <tibble [535 × 2]>
#> 3 河北省 <tibble [747 × 2]>
#> 4 山西省 <tibble [519 × 2]>
#> 5 内蒙古自治区 <tibble [372 × 2]>
#> 6 辽宁省 <tibble [1,092 × 2]>
#> 7 吉林省 <tibble [636 × 2]>
#> 8 黑龙江省 <tibble [500 × 2]>
#> 9 上海市 <tibble [3,536 × 2]>
#> 10 江苏省 <tibble [3,226 × 2]>
#> # … with 25 more rows

nestdf 是个 nested tibble 数据框,下面我们再把这个数据框 unnest:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
nestdf %>% 
unnest(data)
#> # A tibble: 37,625 x 3
#> name code year
#> <chr> <dbl> <dbl>
#> 1 北京市 8 2015
#> 2 北京市 8 2016
#> 3 北京市 8 2017
#> 4 北京市 10 2010
#> 5 北京市 10 2011
#> 6 北京市 10 2013
#> 7 北京市 10 2012
#> 8 北京市 10 2015
#> 9 北京市 10 2014
#> 10 北京市 18 2015
#> # … with 37,615 more rows

输出为 xlsx 文件:

1
2
3
nestdf %>% 
unnest(data) %>%
writexl::write_xlsx("省级分类.xlsx")

琼琼还想把深圳的观测值筛选出来,这就很容易了:

1
2
3
4
# 深圳的
sz <- read_sf('深圳市.json')
szdf <- df %>%
st_intersection(sz)

绘图展示下:

1
2
3
4
5
6
7
8
9
10
11
12
13
ggplot() + 
geom_sf(data = sz, size = 0.5, fill = NA,
color = "white") +
geom_sf_text(data = szdf,
label = "", family = "FontAwesome",
alpha = 0.4, size = 2,
color = "#E31A1C") +
coord_sf(crs = "+proj=laea +lat_0=23 +lon_0=113 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs") +
theme_modern_rc(base_family = cnfont) +
worldtilegrid::theme_enhance_wtg() +
labs(x = "", y = "",
title = "深圳市的上市公司",
subtitle = "TidyFriday Project")

再输出为 xlsx 文件:

1
2
3
4
5
szdf %>% 
select(code, year, name = name.1) %>%
mutate(city = "深圳") %>%
st_drop_geometry() %>%
writexl::write_xlsx("深圳的.xlsx")

这样我们就完成了琼琼的所有要求了。

不过上面这样做是没有意义的。

因为我们可以直接爬取到所有上市公司基本信息的数据,这个数据接口在这里:和讯财经

我之前介绍的 finance 命令包(Stata 命令,建议从知识星球下载源码进行本地安装)中有一个 cnstock2 命令可以用于爬取这个接口(仅限 Mac 用户使用):

1
2
3
4
5
6
7
8
* 下载两市所有上市公司的基本情况数据:
cnstock2
save code.dta, replace

* 下载沪市所有上市公司的基本情况数据:
cnstock2, m(SH)
* 下载深市所有上市公司的基本情况数据:
cnstock2, m(SZ)

code.dta 是这样的:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
library(readstata13)
read.dta13("code.dta") %>%
as_tibble()

#> # A tibble: 3,704 x 12
#> number code name total_stock_num outstanding_sto… outstanding_sto…
#> <int> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 1 6013… 工商银行… 3564. 2696. 14775.
#> 2 2 6005… 贵州茅台… 12.6 12.6 13774.
#> 3 3 6012… 农业银行… 3500. 2941. 10233.
#> 4 4 6013… 中国平安… 183. 108. 8850.
#> 5 5 6018… 中国石油… 1830. 1619. 8630.
#> 6 6 6019… 中国银行… 2944. 2108. 7588.
#> 7 7 6000… 招商银行… 252. 206. 7329.
#> 8 8 6016… 中国人寿… 283. 208. 6495.
#> 9 9 0008… 五粮液… 38.8 38.0 4783.
#> 10 10 6000… 中国石化… 1211. 956. 4529.
#> # … with 3,694 more rows, and 6 more variables:
#> # registered_capital <dbl>, pe_ratio <dbl>, industry <chr>,
#> # concept <chr>, area <chr>, close <dbl>

这里我使用的是 readstata13 这个包读取 dta 文件,这个包的 read.dta13 函数据说是读取 Stata 文件最快速的方式。

Windows 用户就不要运行这个 cnstock2 命令了,不支持。可以直接从知识星球下载附件获取这份数据集。

# R

评论

Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×