Rui - 地图可视化（四）：一些问题

所有企业的分布图

在第三节中，我挑选了频数最高的四类企业/研究院进行可视化。这一次试着将所有企业/研究院的地理分布展现出来。

Code

library(tidyverse)

data <- read.csv("data/地址经纬度数据.csv", encoding = 'UTF-8')
data <- data %>%
  filter(!is.na(longitude), !is.na(latitude))

area_df <- read.csv("越城区边界.csv")

Code

library(leaflet)
library(leafletCN)

Code

m <- leaflet() %>% 
  setView(lng = median(data$longitude), lat = median(data$latitude), zoom = 11) %>%
  amap() %>% 
  addCircleMarkers(
    data = data, 
    lng = ~longitude, 
    lat = ~latitude,
    color = "red",
    stroke = FALSE,
    radius = 4,
    fillOpacity = 0.3
  ) %>%
  addPolygons(
    lng = area_df$longitude,
    lat = area_df$latitude,
    fillColor = "royalblue",
    fillOpacity = 0.2,
    color = "blue",
    weight = 3
  )

m

发现问题

如果将地图缩小可以发现：不少企业被标在了绍兴市越城区以外的地方，甚至标在了云南、广东等地。这明显与事实不符，可能是在使用高德 API 进行地址经纬度转换时导致的。

那么如何判断企业经纬度坐标是否在指定的区域（绍兴市越城区）内呢？可以将其抽象为一个数学问题：判断某点是否在多边形内部。对于这个数学问题，可以用射线法解决：

首先，将多边形的所有顶点按照任意方向进行排序。
选择一个点 P，作为待判断的点。
从点 P 向任意方向发出一条射线，并计算该射线与多边形的每一条边的交点数。
如果交点数为奇数，点 P 在多边形内部；如果交点数为偶数，点 P 在多边形外部。

ChatGPT 也给出了射线法的一个小示例：

Code

# 定义多边形顶点经纬度坐标
poly_coords <- data.frame(
  lat = c(38.894134, 38.895445, 38.895375, 38.894124),
  lng = c(-77.037212, -77.036520, -77.033418, -77.035202)
)

# 定义待判断的点坐标
point_coord <- data.frame(
  lat = 38.894586,
  lng = -77.035177
)

# 判断点是否在多边形内
is_point_in_poly <- function(poly_coords, point_coord) {
  nvert <- nrow(poly_coords)
  testx <- point_coord$lng
  testy <- point_coord$lat
  j <- nvert
  c <- FALSE
  for (i in 1:nvert) {
    if (((poly_coords[i, "lat"] > testy) != (poly_coords[j, "lat"] > testy)) &&
        (testx < (poly_coords[j, "lng"] - poly_coords[i, "lng"]) * (testy - poly_coords[i, "lat"]) /
         (poly_coords[j, "lat"] - poly_coords[i, "lat"]) + poly_coords[i, "lng"])) {
      c <- !c
    }
    j <- i
  }
  return(c)
}

# 判断点是否在多边形内
is_in_poly <- is_point_in_poly(poly_coords, point_coord)

if (is_in_poly) {
  print("Point is inside the polygon.")
} else {
  print("Point is outside the polygon.")
}
## [1] "Point is inside the polygon."

但是俗话说得好：“包治百病”。可以利用现成的包中的某些函数直接实现我们的需求。

我使用的是 sp 包中的 point.in.polygon() 来判断点是否在指定多变形区域内，并且 point.in.polygon() 是向量化的函数。

Code

library(sp)

data$in_polygon <- point.in.polygon(
  point.x = data$longitude, 
  point.y = data$latitude, 
  pol.x = area_df$longitude, 
  pol.y = area_df$latitude
)

data %>% head() %>% knitr::kable()

company_type	name	clean_address	longitude	latitude	in_polygon
共建研究院	浙江大学绍兴研究院	浙江省绍兴市越城区迪荡湖隧道科学园3号楼26-27层	120.6082	30.01323	1
共建研究院	天津大学浙江国际创新设计与智造研究院	浙江省绍兴市越城区洋泾湖科创园1号楼	120.6354	30.05344	1
共建研究院	上海大学绍兴研究院	浙江省绍兴市越城区三江路78号	120.6241	30.10052	1
共建研究院	浙江工业大学绍兴研究院	浙江省绍兴市越城区洋泾湖科创园2号楼	120.6354	30.05344	1
共建研究院	江南大学(绍兴)产业技术研究院	浙江省绍兴市越城区洋江东路19号	120.6075	30.04988	1
共建研究院	中国纺织科学研究院江南分院	浙江省绍兴市越城区双堰路30号	120.6113	30.09832	1

Note

integer array; values are: 0: point is strictly exterior to pol; 1: point is strictly interior to pol; 2: point lies on the relative interior of an edge of pol; 3: point is a vertex of pol.

过滤数据：

Code

data <- data %>%
  filter(in_polygon == 1)

作图：

Code

m <- leaflet() %>% 
  setView(lng = median(data$longitude), lat = median(data$latitude), zoom = 11) %>%
  amap() %>% 
  addCircleMarkers(
    data = data, 
    lng = ~longitude, 
    lat = ~latitude,
    color = "red",
    stroke = FALSE,
    radius = 4,
    fillOpacity = 0.3
  ) %>%
  addPolygons(
    lng = area_df$longitude,
    lat = area_df$latitude,
    fillColor = "royalblue",
    fillOpacity = 0.2,
    color = "blue",
    weight = 3
  )

m