一、 背景技术

1、POI简介

兴趣点(Point of Interest,POI)也可以被称为信息点(Point of Information)，是指地图上具有特定信息的位置点，可以是楼宇，小区，商店，银行，学校等。在现代电子地图中，兴趣点通常都被带有特殊符号的标志标识出来，方便人们找到这些点。POI是基于位置的服务的核心数据，通过POI可以获得丰富的信息，例如用户经常光顾的地点类别，因此很多公司都会根据用户的POI来做数据挖掘，完善用户画像。在电子地图系统中，往往会对POI做完整详细的分类，按照分类可以将POI分为几个类别，例如美食、购物、房地产和旅游景点等。

2、泰森多边形

泰森多边形又叫冯洛诺伊图（Voronoi diagram），得名于Georgy Voronoi，是一组由连接两邻点线段的垂直平分线组成的连续多边形组成。一个泰森多边形内的任一点到构成该多边形的控制点的距离小于到其他多边形控制点的距离。泰森多边形的特性符合基站分布的特点，适合做面状区域的划分。呈现出的多边形可以用来计算基站覆盖面积、半径等，也可以作为进一步分析的基础图层。

二、 详细技术

1、在线地图POI数据获取

城市化的快速发展使城市的每个部分分化出特定的职能，城市中具有特定职能的区域叫做功能区，例如居住区，商业区，工业区等。用户在不同的功能区表现出不同的移动特点，人们在一天的移动过程中也会经过不同的功能区，因此，研究用户在功能区内和功能区之间的移动行为，有利于进一步对用户的移动特征进行分析。研究用户的移动轨迹，结合不同的城市功能区域对人们的移动行为的影响，更进一步探索人们的兴趣区域和移动行为模式。

目前大多数电子地图提供开放的HTTP接口，方便使用者通过这些接口使用各类型的地理数据服务，返回相应的数据格式如XML和JSON供使用者查询。比如百度地图API提供地点检索服务（又名Place API）是一类Web API接口服务。服务提供多种场景的地点（POI）检索功能，包括城市检索、圆形区域检索、矩形区域检索。开发者可通过接口获取地点（POI）基础或详细地理信息。

POI分类

一个请求的示例如下:

http://api.map.baidu.com/place/v2/search?query=黄鹤楼®ion=武汉&output=json&ak=您的ak //GET请求

http://api.map.baidu.com/place/v2/detail?uid=21f7a9ad792605bd9e9d68c9&output=json&scope=2&ak=您的ak //GET请求

获得的一个POI信息如下:

    根据POI类型对地点进行功能区划分，只取名称(name)，类型(type)，地址(address)，行政区划(adname)，位置(location)五个属性来表示。表是过滤字段之后的POI数据的格式。

2、基站区域划分

理想的基站覆盖应该呈蜂窝的形状，并且假设每个基站的覆盖面积都相同，这样可以保证基站的覆盖最大化且无重叠。但在实际中，基站所用的天线往往是全向天线或三向天线，基站的覆盖从平面上看是圆形，这样就会造成基站的覆盖有重叠的部分。在密集城区中，由于用户量多，布站密度大，因此这种基站覆盖相互重叠也更加严重。正是由于这种基站重叠覆盖，会导致用户在两个基站的连接频繁切换。此外，郊区的基站以宏站为主，覆盖半径大，这就使基站网络的拓扑关系变得难以描述。

为了更加准确的描述基站的覆盖范围，采用泰森多边形来表示基站的覆盖面积。将平面划分成以每个基站为中心的Voronoi区域，首先使用Delaunay三角剖分算法构建三角形网，然后遍历每个顶点，计算与之相邻的三角形的外接圆圆心，连接所有外接圆的圆心，就可以得到该顶点的Voronoi区域。在下一小节中，为了计算每个区域内的POI的数量，需要计算每个区域的覆盖范围。由于Voronoi区域是不规则图形，覆盖范围不好量化，本文采用的量化方法为:以该基站作为圆心，圆心到该Voronoi区域最远的顶点的距离作为覆盖半径，所作的圆形作为该区域的覆盖范围。算法4-1是Voronoi区域的形成及覆盖半径计算方法，输入所有基站点的坐标，输出每个基站对应的Vorono区域的外接圆的半径。

Delaunay三角网的生成

（1）构造一个超级三角形，包含所有散点，放入三角形链表。

（2）将点集中的散点依次插入，在三角形链表中找出其外接圆包含插入点的三角形（称为该点的影响三角形），删除影响三角形的公共边，将插入点同影响三角形的全部顶点连接起来，从而完成一个点在Delaunay三角形链表中的插入。

（3）根据优化准则对局部新形成的三角形进行优化。将形成的三角形放入Delaunay三角形链表。

（4）循环执行上述第2步，直到所有散点插入完毕。

关键步骤2如下图所示：

下图展示了某市TD-LTE基站生成的泰森多边形图，其中每一个多边形表示基站的覆盖范围。

3、区域功能区标记

上一节中对基站进行区域划分，本节需要对基站区域进行功能区标记，进而实现对用户轨迹的语义化建模。通过百度地图提供的API接口获取了A市所有的POI数据，其中城区有**条，郊区有**条，共计**条。

本文主要研究城市居民的移动行为，城区是高度发达的城市区域，城市功能健全，城市布局也相对完善。城市中最多的功能区是居住区和办公区，因为大多数城市居民都是“日出而作，日落而息”。本文选取居住区，商业区，办公区，科教区和风景区五个功能区域作为研究的对象。百度地图POI分为两类，分别是一级分类和二级分类，一级分类是POI所属类别，二级分类为细分类别。选取一级分类作为功能区的主要分类依据，对应关系如表所示。

用户的位置是通过基站的位置反映，因此需要将基站区域进行功能区标注。在上一节中，己经将地图按照基站划分为若干的Voronoi区域，并且得到了每个区域的覆盖半径。对于每一个Voronoi区域而言，可能会有多个POI落入该区域内，通过计算落入每个区域内的各POI种类的数量，可以获得每个区域内部不同类别POI所占的比例，根据比例可以判断该区域属于哪一种功能区。对于每一个区域Vor_i，用一个包含不同类别POI的数量的集合来表示，即

基于以上思路，具体的标记步骤如下:

第一步，选定一个基站;

第二步，遍历所有POI，计算POI到基站的距离，若小于基站覆盖半径，则认为此POI属于该基站区域，记录POI的种类并更新POI的数量;若大于基站的覆盖半径，则此POI不属于该基站区域，计算下一个POI与基站的距离。直到遍历完所有POI;

第三步，计算该基站各类别POI的数量，选取数量最多的POI类别作为该基站区域的种类。若存在多个相同数量的种类，计算各种类的质心，选择质心到基站最近的种类作为基站区域种类;

第四步，遍历所有基站，直到计算出所有基站区域的种类。

基站区域功能区分类的算法流程如图所

该市共有TD-LTE小区**个，经过合并相同位置之后的小区**个。标注居住区类基站**个，商业区类基站**个，办公区类基站**个，科教区类基站**个，旅游区类基站**个。标注后的基站区域如图所示。通过对功能区进行标记，便可根据用户在功能区之间的移动顺序对用户轨迹进行语义化建模。

附python取高德POI数据代码

主代码

from urllib.parse import quote

import json

import os

from transCoordinateSystem import gcj02_to_wgs84, gcj02_to_bd09

import area_boundary as  area_boundary

import city_grid as city_grid

import time

import collections

import pandas as pd

from requests.adapters import HTTPAdapter

import requests

#from shp import trans_point_to_shp

#################################################需要修改###########################################################

## TODO 1.划分的网格距离，0.02-0.05最佳，建议如果是数量比较多的用0.01或0.02，如餐厅，企业。数据量少的用0.05或者更大，如大学

pology_split_distance = 0.5

## TODO 2. 城市编码，参见高德城市编码表，注意需要用adcode列的编码

city_code = '150100'

## TODO 3. POI类型编码，类型名或者编码都行，具体参见《高德地图POI分类编码表.xlsx》

typs = ['加油站']  # ['企业', '公园', '广场', '风景名胜', '小学']

## TODO 4. 高德开放平台密钥

gaode_key = ['高德秘钥1', '高德秘钥2']

# TODO 5.输出数据坐标系,1为高德GCJ20坐标系，2WGS84坐标系，3百度BD09坐标系

coord = 2

############################################以下不需要动#######################################################################

poi_pology_search_url = 'https://restapi.amap.com/v3/place/polygon'

buffer_keys = collections.deque(maxlen=len(gaode_key))

def init_queen():

    for i in range(len(gaode_key)):

        buffer_keys.append(gaode_key[i])

    print('当前可供使用的高德密钥：', buffer_keys)

# 根据城市名称和分类关键字获取poi数据

def getpois(grids, keywords):

    if buffer_keys.maxlen == 0:

        print('密钥已经用尽，程序退出！！！！！！！！！！！！！！！')

        exit(0)

    amap_key = buffer_keys[0]  # 总是获取队列中的第一个密钥

    i = 1

    poilist = []

    while True:  # 使用while循环不断分页获取数据

        result = getpoi_page(grids, keywords, i, amap_key)

        print(result)

        print("当前爬取结果:", result)

        if result != None:

            result = json.loads(result)  # 将字符串转换为json

            try:

                if result['count'] == '0':

                    break

            except Exception as e:

                print('出现异常：', e)

            if result['infocode'] == '10001' or result['infocode'] == '10003':

                print(result)

                print('无效的密钥！！！！！！！！！！！！！，重新切换密钥进行爬取')

                buffer_keys.remove(buffer_keys[0])

                try:

                    amap_key = buffer_keys[0]  # 总是获取队列中的第一个密钥

                except Exception as e:

                    print('密钥已经用尽，程序退出...')

                    exit(0)

                result = getpoi_page(grids, keywords, i, amap_key)

                result = json.loads(result)

            hand(poilist, result)

        i = i + 1

    return poilist

# 数据写入csv文件中

def write_to_csv(poilist, citycode, classfield, coord):

    data_csv = {}

    lons, lats, names, addresss, pnames, citynames, business_areas, types, typecodes, ids, type_1s, type_2s, type_3s, type_4s = [], [], [], [], [], [], [], [], [], [], [], [], [], []

    if len(poilist) == 0:

        print("处理完成，当前citycode:" + str(citycode), ", classfield为：", str(classfield) + "，数据为空，，，结束.......")

        return None, None

    for i in range(len(poilist)):

        location = poilist[i].get('location')

        name = poilist[i].get('name')

        address = poilist[i].get('address')

        pname = poilist[i].get('pname')

        cityname = poilist[i].get('cityname')

        business_area = poilist[i].get('business_area')

        type = poilist[i].get('type')

        typecode = poilist[i].get('typecode')

        lng = str(location).split(",")[0]

        lat = str(location).split(",")[1]

        id = poilist[i].get('id')

        if (coord == 2):

            result = gcj02_to_wgs84(float(lng), float(lat))

            lng = result[0]

            lat = result[1]

        if (coord == 3):

            result = gcj02_to_bd09(float(lng), float(lat))

            lng = result[0]

            lat = result[1]

        type_1, type_2, type_3, type_4 = '','','',''

        if str(type) != None and str(type) != '':

            type_strs = type.split(';')

            for i in range(len(type_strs)):

                ty = type_strs[i]

                if i == 0:

                    type_1 = ty

                elif i == 1:

                    type_2 = ty

                elif i == 2:

                    type_3 = ty

                elif i == 3:

                    type_4 = ty

        lons.append(lng)

        lats.append(lat)

        names.append(name)

        addresss.append(address)

        pnames.append(pname)

        citynames.append(cityname)

        if business_area == []:

            business_area = ''

        business_areas.append(business_area)

        types.append(type)

        typecodes.append(typecode)

        ids.append(id)

        type_1s.append(type_1)

        type_2s.append(type_2)

        type_3s.append(type_3)

        type_4s.append(type_4)

    data_csv['lon'], data_csv['lat'], data_csv['name'], data_csv['address'], data_csv['pname'],

    data_csv['cityname'], data_csv['business_area'], data_csv['type'], data_csv['typecode'], data_csv['id'], data_csv[

        'type1'], data_csv['type2'], data_csv['type3'], data_csv['type4'] =

        lons, lats, names, addresss, pnames, citynames, business_areas, types, typecodes, ids, type_1s, type_2s, type_3s, type_4s

    df = pd.DataFrame(data_csv)

    folder_name = 'poi-' + citycode + "-" + classfield

    folder_name_full = 'data' + os.sep + folder_name + os.sep

    if os.path.exists(folder_name_full) is False:

        os.makedirs(folder_name_full)

    file_name = 'poi-' + citycode + "-" + classfield + ".csv"

    file_path = folder_name_full + file_name

    df.to_csv(file_path, index=False, encoding='utf_8_sig')

    print('写入成功')

    return folder_name_full, file_name

# 将返回的poi数据装入集合返回

def hand(poilist, result):

    # result = json.loads(result)  # 将字符串转换为json

    pois = result['pois']

    for i in range(len(pois)):

        poilist.append(pois[i])

# 单页获取pois

def getpoi_page(grids, types, page, key):

    polygon = str(grids[0]) + "," + str(grids[1]) + "|" + str(grids[2]) + "," + str(grids[3])

    req_url = poi_pology_search_url + "?key=" + key + '&extensions=all&types=' + quote(

        types) + '&polygon=' + polygon + '&offset=25' + '&page=' + str(

        page) + '&output=json'

    print('请求url：', req_url)

    s = requests.Session()

    s.mount('http://', HTTPAdapter(max_retries=5))

    s.mount('https://', HTTPAdapter(max_retries=5))

    try:

        data = s.get(req_url, timeout=5)

        return data.text

    except requests.exceptions.RequestException as e:

        data = s.get(req_url, timeout=5)

        return data.text

    return None

def get_drids(min_lng, max_lat, max_lng, min_lat, keyword, key, pology_split_distance, all_grids):

    grids_lib = city_grid.generate_grids(min_lng, max_lat, max_lng, min_lat, pology_split_distance)

    print('划分后的网格数：', len(grids_lib))

    print(grids_lib)

    # 3. 根据生成的网格爬取数据，验证网格大小是否合适，如果不合适的话，需要继续切分网格

    for grid in grids_lib:

        one_pology_data = getpoi_page(grid, keyword, 1, key)

        data = json.loads(one_pology_data)

        print(data)

        while int(data['count']) > 890:

            get_drids(grid[0], grid[1], grid[2], grid[3], keyword, key, pology_split_distance / 2, all_grids)

        all_grids.append(grid)

    return all_grids

def get_data(city, keyword, coord):

    # 1. 获取城市边界的最大、最小经纬度

    amap_key = buffer_keys[0]  # 总是获取队列中的第一个密钥

    max_lng, min_lng, max_lat, min_lat = area_boundary.getlnglat(city, amap_key)

    print('当前城市：', city, "max_lng, min_lng, max_lat, min_lat：", max_lng, min_lng, max_lat, min_lat)

    # 2. 生成网格切片格式：

    grids_lib = city_grid.generate_grids(min_lng, max_lat, max_lng, min_lat, pology_split_distance)

    print('划分后的网格数：', len(grids_lib))

    print(grids_lib)

    all_data = []

    begin_time = time.time()

    print('==========================正式开始爬取啦！！！！！！！！！！！================================')

    for grid in grids_lib:

        # grid格式：[112.23, 23.23, 112.24, 23.22]

        one_pology_data = getpois(grid, keyword)

        print('===================================当前矩形范围：', grid, '总共：',

              str(len(one_pology_data)) + "条数据.............................")

        all_data.extend(one_pology_data)

    end_time = time.time()

    print('全部：', str(len(grids_lib)) + '个矩形范围', '总的', str(len(all_data)), '条数据, 耗时：', str(end_time - begin_time),

          '正在写入CSV文件中')

    file_folder, file_name = write_to_csv(all_data, city, keyword, coord)

    # 写入shp

    #if file_folder is not None:

        #trans_point_to_shp(file_folder, file_name, 0, 1, pology_split_distance, keyword)

if __name__ == '__main__':

    # 初始化密钥队列

    init_queen()

    for type in typs:

        get_data(city_code, type, coord)

transCoordinateSystem代码

#经纬度坐标（WGS84）国测局火星坐标（GCJ02）百度坐标（BD09）相互转换

'''

目前的的坐标体系分类：

一是GPS坐标，也即WGS-84坐标是一个国际的标准，一般卫星导航，原始的GPS设备中的数据都是采用这一坐标系。国外的Google地图、OSM等采用的都是这一坐标。

二是国测局坐标，国测局坐标GCJ-02坐标也叫火星坐标，是国家测绘局为了国家安全在原始坐标的基础上进行偏移得到的坐标，基本国内的电子地图、导航设备都是采用的这一坐标系，如：高德、搜搜、51地图MapABC地图，谷歌中国地图也是。

三是百度坐标，百度坐标BD-09坐标是百度公司出于商业保护在国测局坐标基础上进行的二次加密。

'''

#! /usr/bin/env python3

# -*- coding:utf-8 -*-

__author = "Brady_Hu"

import math

x_pi = 3.14159265358979324 * 3000.0 / 180.0

pi = 3.1415926535897932384626  # π

a = 6378245.0  # 长半轴

ee = 0.00669342162296594323  # 扁率

def gcj02_to_bd09(lng, lat):

    """

    火星坐标系(GCJ-02)转百度坐标系(BD-09)

    谷歌、高德——>百度

    :param lng:火星坐标经度

    :param lat:火星坐标纬度

    :return:

    """

    z = math.sqrt(lng * lng + lat * lat) + 0.00002 * math.sin(lat * x_pi)

    theta = math.atan2(lat, lng) + 0.000003 * math.cos(lng * x_pi)

    bd_lng = z * math.cos(theta) + 0.0065

    bd_lat = z * math.sin(theta) + 0.006

    return [bd_lng, bd_lat]

def bd09_to_gcj02(bd_lon, bd_lat):

    """

    百度坐标系(BD-09)转火星坐标系(GCJ-02)

    百度——>谷歌、高德

    :param bd_lat:百度坐标纬度

    :param bd_lon:百度坐标经度

    :return:转换后的坐标列表形式

    """

    x = bd_lon - 0.0065

    y = bd_lat - 0.006

    z = math.sqrt(x * x + y * y) - 0.00002 * math.sin(y * x_pi)

    theta = math.atan2(y, x) - 0.000003 * math.cos(x * x_pi)

    gg_lng = z * math.cos(theta)

    gg_lat = z * math.sin(theta)

    return [gg_lng, gg_lat]

def wgs84_to_gcj02(lng, lat):

    """

    WGS84转GCJ02(火星坐标系)

    :param lng:WGS84坐标系的经度

    :param lat:WGS84坐标系的纬度

    :return:

    """

    if out_of_china(lng, lat):  # 判断是否在国内

        return lng, lat

    dlat = _transformlat(lng - 105.0, lat - 35.0)

    dlng = _transformlng(lng - 105.0, lat - 35.0)

    radlat = lat / 180.0 * pi

    magic = math.sin(radlat)

    magic = 1 - ee * magic * magic

    sqrtmagic = math.sqrt(magic)

    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)

    dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi)

    mglat = lat + dlat

    mglng = lng + dlng

    return [mglng, mglat]

def gcj02_to_wgs84(lng, lat):

    """

    GCJ02(火星坐标系)转GPS84

    :param lng:火星坐标系的经度

    :param lat:火星坐标系纬度

    :return:

    """

    if out_of_china(lng, lat):

        return lng, lat

    dlat = _transformlat(lng - 105.0, lat - 35.0)

    dlng = _transformlng(lng - 105.0, lat - 35.0)

    radlat = lat / 180.0 * pi

    magic = math.sin(radlat)

    magic = 1 - ee * magic * magic

    sqrtmagic = math.sqrt(magic)

    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)

    dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi)

    mglat = lat + dlat

    mglng = lng + dlng

    return [lng * 2 - mglng, lat * 2 - mglat]

def bd09_to_wgs84(bd_lon, bd_lat):

    lon, lat = bd09_to_gcj02(bd_lon, bd_lat)

    return gcj02_to_wgs84(lon, lat)

def wgs84_to_bd09(lon, lat):

    lon, lat = wgs84_to_gcj02(lon, lat)

    return gcj02_to_bd09(lon, lat)

def _transformlat(lng, lat):

    ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat +

          0.1 * lng * lat + 0.2 * math.sqrt(math.fabs(lng))

    ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *

            math.sin(2.0 * lng * pi)) * 2.0 / 3.0

    ret += (20.0 * math.sin(lat * pi) + 40.0 *

            math.sin(lat / 3.0 * pi)) * 2.0 / 3.0

    ret += (160.0 * math.sin(lat / 12.0 * pi) + 320 *

            math.sin(lat * pi / 30.0)) * 2.0 / 3.0

    return ret

def _transformlng(lng, lat):

    ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng +

          0.1 * lng * lat + 0.1 * math.sqrt(math.fabs(lng))

    ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *

            math.sin(2.0 * lng * pi)) * 2.0 / 3.0

    ret += (20.0 * math.sin(lng * pi) + 40.0 *

            math.sin(lng / 3.0 * pi)) * 2.0 / 3.0

    ret += (150.0 * math.sin(lng / 12.0 * pi) + 300.0 *

            math.sin(lng / 30.0 * pi)) * 2.0 / 3.0

    return ret

def out_of_china(lng, lat):

    """

    判断是否在国内，不在国内不做偏移

    :param lng:

    :param lat:

    :return:

    """

    return not (lng > 73.66 and lng < 135.05 and lat > 3.86 and lat < 53.55)

if __name__ == '__main__':

    lng = 113.224367

    lat = 25.69346

    #result1 = gcj02_to_bd09(lng, lat)

    #result2 = bd09_to_gcj02(lng, lat)

    #result3 = wgs84_to_gcj02(lng, lat)

    #result4 = gcj02_to_wgs84(lng, lat)

    result5 = bd09_to_wgs84(lng, lat)

    #result6 = wgs84_to_bd09(lng, lat)

print(result5)

area_boundary代码

# -*- coding:utf-8 -*-

# 第一行必须有，否则报中文字符非ascii码错误

import urllib.request

from urllib.parse import quote

import json

url = 'http://restapi.amap.com/v3/config/district?'

def getlnglat(address, key):

    uri = url + 'keywords=' + quote(address) + '&key=' + key + '&subdistrict=1' + '&extensions=all'

    print(uri)

    # 访问链接后，api会回传给一个json格式的数据

    temp = urllib.request.urlopen(uri)

    temp = json.loads(temp.read())

    # polyline是坐标，name是区域的名字

    Data = temp["districts"][0]['polyline']

    lngs = []

    lats = []

    points = []

    for line in str(Data).split(";"):

        if len(line.split("|")) > 1:

            for uu in line.split("|"):

                if float(uu.split(",")[0]) != None:

                    lngs.append(float(uu.split(",")[0]))

                    lats.append(float(uu.split(",")[1]))

                    points.append([float(uu.split(",")[0]), float(uu.split(",")[1])])

        else:

            if float(line.split(",")[0]) != None:

                lngs.append(float(line.split(",")[0]))

                lats.append(float(line.split(",")[1]))

                points.append([float(line.split(",")[0]), float(line.split(",")[1])])

    print(points)

    print(max(lngs), min(lngs), max(lats), min(lats))

    return max(lngs), min(lngs), max(lats), min(lats)

'''

num = 0

#ad = getSubName(addr_name)  # 得到福州下属区域的城市代码

add = getlnglat(addr_name)  # 得到福州整个的边界数据

while num < len(ad):

    add = pd.concat([add, getlnglat(ad[num].encode("utf-8"))])  # 得到福州下属的全部区域的边界数据

    num += 1

add.to_csv('{0}.csv'.format(addr_name), encoding='gbk')

'''

city_grid代码

import numpy as np

def generate_grids(start_long,start_lat,end_long,end_lat,resolution):

    """

    根据起始的经纬度和分辨率，生成需要需要的网格.

    方向为左上，右下，所以resolution应为 负数，否则未空

    :param start_long:

    :param start_lat:

    :param end_long:

    :param end_lat:

    :param resolution:

    :return:

    """

    assert start_long < end_long,'需要从左上到右下设置经度，start的经度应小于end的经度'

    assert start_lat > end_lat,'需要从左上到右下设置纬度，start的纬度应大于end的纬度'

    assert resolution>0,'resolution应大于0'

    #assert resolution2>0,'resolution应大于0'

    grids_lib=[]  #空列表

    longs = np.arange(start_long,end_long,resolution)

    if longs[-1] != end_long:

        longs = np.append(longs,end_long)

    lats = np.arange(start_lat,end_lat,-resolution)

    if lats[-1] != end_lat:

        lats = np.append(lats,end_lat)

    for i in range(len(longs)-1):

        for j in range(len(lats)-1):

            grids_lib.append([round(float(longs[i]),6),round(float(lats[j]),6),round(float(longs[i+1]),6),round(float(lats[j+1]),6)])

            #yield [round(float(longs[i]),6),round(float(lats[j]),6),round(float(longs[i+1]),6),round(float(lats[j+1]),6)]

    return grids_lib

"""

将数据存储

grids_lib = generate_grids(114.237373, 30.687372, 114.460680, 30.419513,0.05)

print('划分后的网格数：', len(grids_lib))

#print(grids_lib)

import csv

import codecs

def data_write_csv(file_name, datas):#file_name为写入CSV文件的路径，datas为要写入数据列表

    file_csv = codecs.open(file_name,'w+','utf-8')#追加

    writer = csv.writer(file_csv, delimiter=' ', quotechar=' ', quoting=csv.QUOTE_MINIMAL)

    for data in datas:

        writer.writerow(data)

    print("保存文件成功，处理结束")

data_write_csv('D:/11.csv',grids_lib)

"""