用编程赋能工作系列——地理围栏的基本操作运算

数据小磨坊

发布于 2020-03-06 14:32:12

2.7K1

发布于 2020-03-06 14:32:12

文章被收录于专栏：数据小魔方数据小魔方

接着上一篇的地图系列相关知识，本篇给大家介绍一种局部空间分析的地理围栏运算，具体场景主要用在分析局部的商圈、商场、街道、步行街内部相关变量方面。

这些区域通常没有标准的行政界线，但是在必要的场合，你又非得在地图上将其边界展示出来，并且判断出那些点是在围栏内部，那些点是在围栏外部。

如下图所示，通过前期调研，假如你已经确认了目标分析区域如图中不规则多边形所示，通过地图围栏围栏可以拿到边界经纬度信息，然后需要甄别出待分析的原始数据中，那些点是在目标分析区域内部，并且单独摘出来进行更加细致的分析。

以上过程存在两个难点，目标区域的边界信息如何获取？有了边界信息我如何对自己的原始数据中的点击进行点归属判断？以下内容就是要重点解决这个问题。

如何获取围栏边界信息？

第一个问题很好解决，只要你自己心里对目标分析区域范围有个大概认知，当前很多在线地图商都提供空间围栏的自助圈选服务，仅需手工操作即可。

以下是高德地图提供的自助围栏服务地址：

https://lbs.amap.com/api/javascript-api/example/overlayers/polygon-draw-and-edit

假如你要分析的目标商圈是王府井，我仅需通过地图平台大概知道王府井的具体位置、大概轮廓，就可以通过高德围栏功能进行围栏信息的获取。

（你可以通过https://lbs.amap.com/console/show/picker获取任意一个点的经纬度，对于王府井地区，可以大致取一个中心点以及三个以上的点组成的简要轮廓）

将这些点、轮廓按照下图红色框内所示填入，并点击右上角的运行按钮，页面即可大致定位到王府井，然后你可以在底部编辑器点击编辑按钮，随意拖拉拽页面上的节点编辑器，直至调整到符合心意位置，并点击退出编辑。

退出编辑后，在浏览器右下角的命令行界面运行如下代码：

console.log('当前围栏是：');
        console.log(path);

即可获得一串目标围栏的经纬度字符串。

如何处理围栏并进行点归属判断？

这问题是主要操作难点，涉及到空间数据操纵，以下仍然是两个工具分别讲解：

R语言中的处理方案：

# 将围栏数据改造成R语言中sf包可识别的形式
# 因为原始围栏是一次将经维度按顺序组合并一次拼接起来的，所以需要使用简单的
# 程序转换为含经度、维度的数据框
library('sf')
library('ggplot2')
library('magrittr')
library("leaflet")
library('stringr')

areaFence = '116.412186,39.910991,116.412351,39.910848,116.412639,39.910781,116.415811,39.910852,116.418102,39.910699,116.418235,39.910615,116.418305,39.910447,116.418331,39.908238,116.418252,39.908026,116.418084,39.907943,116.411998,39.907690,116.410818,39.907655,116.410649,39.907764,116.410575,39.907951,116.410543,39.911500,116.410474,39.911654,116.410364,39.911747,116.408000,39.911651,116.407846,39.911743,116.407758,39.911908,116.407709,39.915676,116.407802,39.915830,116.408023,39.915918,116.410122,39.915942,116.410287,39.916017,116.410341,39.916205,116.410104,39.924552,116.410142,39.924653,116.410293,39.924703,116.411451,39.924741,116.411586,39.924687,116.411657,39.924571,116.411881,39.919409,116.412098,39.911179,116.412186,39.910991'
my_fun <- function(data){
  range_data = str_split(data,',')[[1]] %>% as.numeric()
  list_a = c()
  list_b = c()
  for (i in 1:length(range_data)){
    if(i%%2 != 0){
      list_a = c(list_a,range_data[i])
    } else {
      list_b = c(list_b,range_data[i])
    }
  }
  result_data = paste0(list_a,' ',list_b)  %>% paste0(collapse = ",")  %>% sprintf('POLYGON ((%s))',.) 
  return(result_data)
}

# 这里将已经转换为数据框的围栏经纬度信息转换为sf模式的多边形对象
ploygon_data <- my_fun(areaFence)
mapdata <- ploygon_data  %>% 
  st_as_sfc() %>%  
  st_cast("POLYGON")
  
 # 围栏基本信息计算
map_data <- st_sf(geom = mapdata,crs = 4326)  # 添加投影
area = map_data %>% st_area() %>% units::set_units( km^2) #获取面积信息
center = st_centroid(mapdata)[[1]] #获取多边形空间中心点信息
bbox = st_bbox(mapdata)  #获取多边形四至信息
radius = st_distance(center,st_point(c(bbox[1],bbox[2])))
#  #获取多边形中心点和外围点半径（我取了一个左下边界点）

现在打印一下我们获取的围栏在地图上的样子。

#打印围栏边界
nc.sp <- as(mapdata,"Spatial") 
leaflet(nc.sp) %>% 
  setView(center[1],center[2], zoom = 14) %>% 
  addTiles('http://webrd02.is.autonavi.com/appmaptile?lang=zh_cn&size=1&scale=1&style=8&x={x}&y={y}&z={z}',tileOptions(tileSize=256, minZoom=9, maxZoom=17),group="高德地图")  %>%
  addPolygons(fillColor = 'blue',weight=1)

围栏有了，接下来伪造一份分析数据，这份数据中的点围绕以上围栏区域的中心和半径随机分布（具体半径会更大）。

# 以围栏中心为圆点，生成特定半径内模拟随机数
get_random_pos <- function(center_x,center_y,radius){
  u = runif(1000, 0, 1) *radius
  v = runif(1000, 0, 1) *2*pi
  x = center_x + u*cos(v)
  y = center_y + u*sin(v)
  return(data.frame(x,y))
}
point_data <- get_random_pos(center_x = 116.412,center_y = 39.91351,radius  = 0.008)
# 打印随机数
 leaflet(nc.sp) %>% 
  setView(center[1],center[2], zoom = 14) %>% 
  addTiles('http://webrd02.is.autonavi.com/appmaptile?lang=zh_cn&size=1&scale=1&style=8&x={x}&y={y}&z={z}',tileOptions(tileSize=256, minZoom=9, maxZoom=17),group="高德地图")  %>%
  addPolygons(fillColor = 'blue',weight=1)  %>%
  addCircleMarkers(
    data = point_data,
    lng = ~x,
    lat = ~y,
    #stroke = FALSE,
    radius = 2,
    fillColor = 'red'
  )

利用sf包中的点归属判别函数st_contains，把随机点都打上一个是否在围栏内部的标签，方便之后在呈现层进行区别填色。

# 判别随机数中那些点落在围栏内部
Point_judge <- function(x,y){
  point_data_a = y
  point_data_b = y  %>% st_as_sf(coords = c("x", "y")) %>%  st_set_crs(4326) %>% st_cast("POINT")
  contain_label =  st_contains(x,point_data_b$geometry) %>% unlist()
  # 核心是这局判别逻辑
  for(i in 1:nrow(y)){
    point_data_a$type[i] = ifelse(i %in% contain_label,'in','out')
  }
  return(point_data_a)
}
points_data <- Point_judge(map_data,point_data)

# 打印可以区分围栏内部&外部的点
cpal <- colorFactor(c("#fb832d",'#098154'), domain = unique(points_data$type)) 
leaflet(nc.sp) %>% 
  setView(center[1],center[2], zoom = 14) %>% 
  addTiles('http://webrd02.is.autonavi.com/appmaptile?lang=zh_cn&size=1&scale=1&style=8&x={x}&y={y}&z={z}',tileOptions(tileSize=256, minZoom=9, maxZoom=17),group="高德地图")  %>%
  addPolygons(fillColor = 'blue',weight=1)  %>%
  addCircleMarkers(
    data = points_data ,
    lng = ~x,
    lat = ~y,
    stroke = FALSE,
    fillColor = ~cpal(type),
    fillOpacity = 0.4,
    radius = 5
  )

可以看到目前属于所圈定的王府井区域内部的点已经被标成了橙黄色，区外的点统一标成了绿色，以上就达到了我们的总体目标，如有后续分析，就可以针对区内的点进行更为深入和细致的分析，比如下钻到非常细微的网格颗粒度进行空间聚类、识别高密度区域等。

Python中的处理方案：

from shapely.geometry import Polygon
from shapely.geometry import Point
from numpy import asarray
import folium
import geopandas as gpd
import pandas as pd
import numpy as np

# 将围栏数据改造成Python中shaplely包可识别的形式
areaFence = '116.412186,39.910991,116.412351,39.910848,116.412639,39.910781,116.415811,39.910852,116.418102,39.910699,116.418235,39.910615,116.418305,39.910447,116.418331,39.908238,116.418252,39.908026,116.418084,39.907943,116.411998,39.907690,116.410818,39.907655,116.410649,39.907764,116.410575,39.907951,116.410543,39.911500,116.410474,39.911654,116.410364,39.911747,116.408000,39.911651,116.407846,39.911743,116.407758,39.911908,116.407709,39.915676,116.407802,39.915830,116.408023,39.915918,116.410122,39.915942,116.410287,39.916017,116.410341,39.916205,116.410104,39.924552,116.410142,39.924653,116.410293,39.924703,116.411451,39.924741,116.411586,39.924687,116.411657,39.924571,116.411881,39.919409,116.412098,39.911179,116.412186,39.910991'
def my_fun(areaFence):
    range_data = [float(i) for  i in  areaFence.split(',')]
    list_a = []
    list_b = []
    for i in range(1,len(range_data)+1):
        if (i%2) != 0:
            list_a.append(range_data[i-1])
        else:
            list_b.append(range_data[i-1])
    return np.vstack([list_a,list_b]).T
my_result = my_fun(areaFence)

# 获取围栏基本信息
polygon=Polygon(my_result)
'{:.10f}'.format(polygon.area)  #围栏面积
centre = tuple(asarray(polygon.centroid).tolist()) #围栏中心点
polygon.bounds # 围栏四至边界点
radius = polygon.centroid.distance(Point(polygon.bounds[0],polygon.bounds[1])) #围栏中心与左下边界距离

使用Python中的folium包来进行打印，这个表也是调用的leaflet在线地图。


#围栏打印
Polygon = gpd.GeoDataFrame( crs= {'init': 'epsg:4326'}, geometry=[polygon]) 

m = folium.Map(
    location = (centre[1],centre[0]), 
    zoom_start=15,
    tiles = 'http://webrd02.is.autonavi.com/appmaptile?lang=zh_cn&size=1&scale=1&style=8&x={x}&y={y}&z={z}',
    attr = '&copy; <a href="http://ditu.amap.com/">高德地图</a>'
)
folium.GeoJson(Polygon).add_to(m)
m

通过简单生成一组围绕围栏中心分布的随机点，来制作一份备用样本。

# 以围栏中心为圆点，生成特定半径内模拟随机数
def get_random_pos(center_x,center_y,radius):
  u = np.random.rand(1000) * radius
  v = np.random.rand(1000) *2*np.pi
  x = centre[0] + np.array(u*np.cos(v))
  y = centre[1] + np.array(u*np.sin(v))
  return np.vstack([y,x])
# 打印随机数
m = folium.Map(
    location = (centre[1],centre[0]), 
    zoom_start=15,
    tiles = 'http://webrd02.is.autonavi.com/appmaptile?lang=zh_cn&size=1&scale=1&style=8&x={x}&y={y}&z={z}',
    attr = '&copy; <a href="http://ditu.amap.com/">高德地图</a>'
)
def out_pic():
    tem_data = get_random_pos(center_x = centre[0],center_y = centre[1],radius  = radius + 0.0015)
    geo_data = [(i,j) for i,j in  zip(tem_data[0],tem_data[1])]
    folium.GeoJson(Polygon).add_to(m)
    for e in geo_data:
        folium.CircleMarker(
        [e[0],e[1]],
        #popup=str(e[2])+":"+str(e[3]),
        fill_color='red',
        color='red',
        fill=True,
        number_of_sides=10,
        radius=2.5,
        stroke = True,
        opacity = 0.5,
        fillOpacity = 1000
        ).add_to(m)
    display(m)
    return geo_data
out_data =  out_pic()

使用Python中shapely包（底层也是和R语言中的sf包基于相同的理论基础实现的）提供的点判别函数contains。

# 判别随机数中那些点落在围栏内部
def Point_judge(out_data,Polygon_data):
    temp = list()
    for i,j in out_data:
        if Polygon_data.contains(Point(j,i))[0]:
            temp.append(1)
        else:
            temp.append(0)
    result_data = np.column_stack((out_data,np.array(temp)))
    return result_data
Points_data = Point_judge(out_data,Polygon)


# 打印可以区分围栏内部&外部的点
m = folium.Map(
    location = (centre[1],centre[0]), 
    zoom_start=15,
    tiles = 'http://webrd02.is.autonavi.com/appmaptile?lang=zh_cn&size=1&scale=1&style=8&x={x}&y={y}&z={z}',
    attr = '&copy; <a href="http://ditu.amap.com/">高德地图</a>'
)
def get_random_pos(center_x,center_y,radius):
  u = np.random.rand(1000) * radius
  v = np.random.rand(1000) *2*np.pi
  x = centre[0] + np.array(u*np.cos(v))
  y = centre[1] + np.array(u*np.sin(v))
  return np.vstack([y,x])
def out_pict():
    tem_data = get_random_pos(center_x = centre[0],center_y = centre[1],radius  = radius + 0.0015)
    folium.GeoJson(Polygon).add_to(m)
    for e in Points_data:
        color_t = str(np.where(e[2] == 1,'red', 'blue')),
        folium.CircleMarker(
        [e[0],e[1]],
        #popup=str(e[2])+":"+str(e[3]),
        fill_color =  color_t,
        color = color_t,
        fill=True,
        number_of_sides=10,
        radius=5,
        stroke = True,
        opacity = 0.5,
        fillOpacity = 1
        ).add_to(m)
    display(m)
out_data =  out_pict()