python 读取micaps第14类数据，并用Basemap进行绘图

18393810639

https://blog.csdn.net/weixin_43718675/article/details/89428532
1. 前言
在大气科学领域，经常用到micaps这个软件。该软件下载链接如下。
http://www.micaps.cn/MiFun/main

在基于micaps进行天气分析时，需要在micaps底图上勾勒出一些线条，用来表示槽线或者辐合线等。如下图所示：

                               
登录/注册后可看大图

将包含这些线条的交互图层进行保存后，就会得到micaps第14类数据。
micaps数据类型介绍如下：https://www.bbsmax.com/A/WpdKqqOrJV/

                               
登录/注册后可看大图

2. 需求
由于在micaps上画出的线条进行保存得到的图片不太符合论文要求。因此有时候需要将这些线条的经纬度坐标记住，然后在利用其他绘图软件进行图形绘制。
如第一张图，共计有4条线，其保存的数据为无格式文本数据，用记事本打开后，实际数据如下。可以看出，每条线的组成点的经纬度都在上面。

                               
登录/注册后可看大图

但是这样看优点不规整，因此可以借助python来读取。

3. python进行读取

```
   filename='../../dryline_position/2015_position/201506'
    f=open(filename,mode='r')
    data=f.readlines()  
```
data的数据格式为：

                               
登录/注册后可看大图

注意，这里的每一行都是字符串。
第一行：交互图层保存的路径
第二行：保存时间
第四行：线条数
第五行：第一条线的描述，其中最后一个数字37表示这条线有多少个点组成。
6-15行为这37个点的描述：经度、维度、和一个不太重要的量。

另外，由于数据行，每行只包括4个数据点，因此37个点，需要int(np.ceil(37/4))，即往上取整数量的行数来表示数据。

第16行：第二条线的描述，同样，最后一个数字为第二条线的组成点数。后面依次类推。


4. 具体代码

```
def get_drylines_position_coordinate(filename):
    '''
    读取micaps第14类数据，micaps数据类型介绍如下链接：
    https://www.bbsmax.com/A/WpdKqqOrJV/

    input: 文件路径
    output: micaps底图上叠加的线条对应的经纬度
    '''
    #读取文件
    f=open(filename,mode='r')
    data=f.readlines()  

    #获取该文件中包含几条线
    num_drylines=int(data[3].strip().split()[-1])

    n=4 #前五行为数据说明

    all_position=[]

    #获取每条线的经纬度坐标
    for i in range(num_drylines):

        #获取该线条用多少个点来描述，其中每个点对应一个经度和维度(和一个额外的数字，不重要)
        n_points=int(data[n].strip().split()[-1])

        #由于默认每行仅支持输出个4数据点，因此，会存在多行数据
        n_rows=int(np.ceil(n_points/4))

        #把包含这条线所有点的列连接起来
        point_data=data[n+1]
        for i in range(n_rows-1):
            point_data=point_data+data[n+2+i]

        #去除换行符和将每个数据隔开
        point_data=point_data.strip().split()

        #将数据类型由str转换为float型
        point_data=[float(point_data) for i in range(len(point_data))]

        #前面说到，每个点包含三个数据，其中第一个为经度，第二个为维度
        point_lon=point_data[0::3]
        point_lat=point_data[1::3]

        #将该条线的经纬度坐标记录下来
        coordinate=[point_lon,point_lat]
        all_position.append(coordinate)

        #每条线差两列
        n=n+n_rows+2

    return all_position ,data
```
下面进行测试：

```
filename='../../dryline_position/2015_position/201506'
coor_data,data=get_drylines_position_coordinate(filename)
```
其中的coor_data为：

                               
登录/注册后可看大图

每条线的经纬度用一个列表来表示。我们点开第一个列表，
发现里面还是一个二维列表。这表示第一条线的组成点的经纬度。
后面依次类推。


后续添加

1. 基于python的basemap库在地图底图上画出上述四条线

```
import numpy as np
import pandas as pd
import os
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
```

```
def plot_map(data,lat_min=35,lat_max=55,lon_min=115,lon_max=135):

    figure=plt.figure(figsize=(16,10))

    #设置投影方式和经纬度范围
    m=Basemap(projection='cyl',llcrnrlat=lat_min,llcrnrlon=lon_min,
          urcrnrlat=lat_max,urcrnrlon=lon_max,resolution='l')

    #画出海岸线和国境线
    m.drawcoastlines()
    m.drawcountries(linewidth=1)   

    #设置坐标刻度
    x_grid=np.arange(lon_min,lon_max+1,5)
    y_grid=np.arange(lat_min,lat_max+1,5)
    m.drawparallels(y_grid)
    m.drawmeridians(x_grid)
    plt.xticks(x_grid,x_grid,fontsize=16)
    plt.yticks(y_grid,y_grid,fontsize=16)
    plt.xlabel('longitude:°E',fontsize=20)
    plt.ylabel('latitude:°N',fontsize=20)

    #读取省界数据。下载省界数据链接为：https://blog.csdn.net/weixin_36677127/article/details/83314583
    m.readshapefile('../../gadm36_CHN_shp/gadm36_CHN_1', 'states',
                    drawbounds = True,linewidth=1)

    #画出每条线
    num_lines=len(coor_data)
    for i in range(num_lines):
        m.plot(data[0],data[1])

    return None
```
下载省界数据链接为：https://blog.csdn.net/weixin_36677127/article/details/83314583
测试一波
```
a=plot_map(coor_data)
```
结果如下图。这种图是不是很学术，哈哈哈哈

                               
登录/注册后可看大图

欢迎交流！

附上测试数据：

tulalang

楼主思路真的太棒了，不过图片没上传成功？我这边看不见图片

18393810639

tulalang 发表于 2019-4-21 10:53
楼主思路真的太棒了，不过图片没上传成功？我这边看不见图片

应该上传成功了，我自己都能看到。如果没有看到图片，你可以点击刚开始的那一个链接，这是我在CSDN上写的。那个编辑器好用点。

tulalang

真的很棒，谢谢楼主分享！