【Basemap 绘制二维网格图】

Basemap 绘制二维网格图

需求：根据四至经纬度绘制规定长的网格图

"""
网格的脚本
"""
from matplotlib import pyplot as plt
from typing import List, Tuple
from mpl_toolkits.basemap import Basemap
import numpy as np
import math

class Error(Exception):
    pass


class ParamInvalid(Error):
    pass


def get_step_array(start: float, end: float, step: float = 500):
    tmp = start
    _arr = []
    while tmp <= end:
        _arr.append(round(tmp, 6))
        tmp += step
    return _arr


def draw_grid(base_map, rect: List[Tuple[float, float]], grid_len: float):
    """
    画网格

    :param base_map: [Basemap] 地图画布
    :param rect: List[Tuple[float, float]] 四至经纬度列表[(lng, lat)]
    :param grid_len: [int] 网格边长
    :return:
    """
    if len(rect) != 4:
        raise ParamInvalid()

    # 获取四至点
    left_top_point, left_down_point, right_top_point, right_down_point = rect
    print(right_top_point)

    # 获取四至点的经纬度
    left_top_lng, left_top_lat = left_top_point
    left_down_lng, left_down_lat = left_down_point
    right_top_lng, right_top_lat = right_top_point
    right_down_lng, right_down_lat = right_down_point

    min_lng = left_top_lng
    max_lng = right_top_lng
    min_lat = left_down_lat
    max_lat = right_top_lat
    print(min_lat, max_lat)
    lons = get_step_array(min_lng, max_lng, step=grid_len)
    lats = get_step_array(min_lat, max_lat, step=grid_len)

    print("经度分割")
    print(f"{lons}")
    print("-------------------")
    print("纬度分割")
    print(lats)
    print(lons)

    return lons, lats



def generate_map(llcrnrlat, llcrnrlon, urcrnrlon, urcrnrlat) -> Basemap:
    """
    生成地图画布

    :return:
    """
    _m = Basemap(llcrnrlat=llcrnrlat, llcrnrlon=llcrnrlon, urcrnrlat=urcrnrlat, urcrnrlon=urcrnrlon)
    return _m


if __name__ == '__main__':
    # 四至坐标
    _rect = [
        # 左上
        (116.14, 39.8),
        # 左下
        (116.14, 39.441),
        # 右上
        (116.785, 39.8),
        # 右下
        (116.785,39.441),
    ]
    #

   
    m = Basemap(llcrnrlon=116.14, llcrnrlat=39.441, urcrnrlon=116.785, urcrnrlat=39.8, lat_1=33,
                 lat_2=45, lon_0=100)
    m.readshapefile(shapefile=r"D:\..\..\..\shapefiles\大兴2d\Export_Output_2",name="Export_Output_2", drawbounds=True)
    # 读取shp文件 我这里用的是2d的，3d的没有试过不知道能不能行

    # 绘制世界地图
    # m.drawcoastlines() #绘制海岸线
    # m.drawcountries()#绘制国界线
    # m.drawmapboundary()#绘制海洋
    # 填充陆地、胡泊、海洋的颜色
    # m.fillcontinents(color='darkolivegreen',  # 陆地颜色
    #                  lake_color='aqua',  # 湖泊颜色#alpha=0.4
    #                   )
    #
    # m.drawmapboundary(fill_color='aqua')  # 填充海洋



    # 区域范围
    bound = {}
    bound['min_lon'] = 116.14 # min_lon
    bound['min_lat'] = 39.441 # min_lat
    bound['max_lon'] = 116.785  # max_lon
    bound['max_lat'] = 39.8  # max_la


    lons, lats =[],[]
    grid_length = 1  # 网格边长是1公里
    mid_lat = (bound['min_lat'] + bound['max_lat']) / 2  # 中间纬度

    #todo 需要考虑不同经纬度所在位置的影响，暂时忽略

    # 如果地球是个规则的球体，其周长是40076千米，那么1个360度的经度圈周长也是40076千米（和赤道圈一样），推断出1度经度的距离是40076 / 360 = 111千米；而1纬度的距离和它所在的纬度带有关。
    #1度就111km
    #具体差距需要考虑不同经纬度所在位置的影响，暂时忽略
    # 1公里约等于0.01度
    lons, lats = draw_grid(base_map=m, rect=_rect, grid_len=0.01)
    m.drawmeridians(lons, linewidth=0.5)
    m.drawparallels(lats, linewidth=0.5)
    plt.show()

基础样式图，后续可以根据需求样式具体调整