美文网首页
如何根据站点污染物数据插值绘制出污染物的空间分布?⑨

如何根据站点污染物数据插值绘制出污染物的空间分布?⑨

作者: Aerosols | 来源:发表于2022-03-31 15:52 被阅读0次

嘿嘿,让我们采取倒叙的方式,先看成品,再看制作过程。

jz.png
  1. 准备焦作市各区县的面图文件,具体操作参考宝藏文章,一步一步地就做好了。
    Arcgis(一) 制作全国行政区shp文件(精确到县级)_皖渝的博客-CSDN博客_arcgis制作shp文件

需要用到全国地理信息库:
成果数据1:100万 (webmap.cn)
需要电脑上有arcgis软件。

  1. 安装一些python库
    参考python geopandas库安装 - 知乎 (zhihu.com)

  2. 照葫芦画瓢
    https://mp.weixin.qq.com/s/shzcZWHLaUOeHGIc7E3Vrg

  3. 我的代码

## 导入库,定义函数
import pandas as pd
import numpy as np
from pykrige.ok import OrdinaryKriging
import plotnine
from plotnine import *
import geopandas as gpd
import shapefile
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
import cmaps
from matplotlib.path import Path
from matplotlib.patches import PathPatch

from scipy.interpolate import griddata
#from cartopy.io.shapereader import Reader
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter

def shp2clip(originfig, ax, shpfile):
    sf = shapefile.Reader(shpfile)
    vertices = []
    codes = []
    for shape_rec in sf.shapeRecords():
        pts = shape_rec.shape.points
        prt = list(shape_rec.shape.parts) + [len(pts)]
        for i in range(len(prt) - 1):
            for j in range(prt[i], prt[i + 1]):
                vertices.append((pts[j][0], pts[j][1]))
            codes += [Path.MOVETO]
            codes += [Path.LINETO] * (prt[i + 1] - prt[i] - 2)
            codes += [Path.CLOSEPOLY]
        clip = Path(vertices, codes)
        clip = PathPatch(clip, transform=ax.transData)
    for contour in originfig.collections:
        contour.set_clip_path(clip)
    return contour

## 读数据
df = pd.read_excel(r"D:\jz_21sta_MonthMean.xlsx",sheet_name="jz_18")

pos1 = pd.read_excel(r"D:\所有空气点位经纬度登记表.xlsx",sheet_name="sta-18")
pos2 = pd.read_excel(r"D:\所有空气点位经纬度登记表.xlsx",sheet_name="sta-12")

## 循环绘图
j= -1
need_time = ["2019/1/31","2019/2/28","2019/3/31","2019/4/30","2019/5/31","2019/6/30","2019/7/31",
             "2019/8/31","2019/9/30","2019/10/31","2019/11/30","2019/12/31","2020/1/31",
             "2020/2/29","2020/3/31","2020/4/30","2020/5/31","2020/6/30","2020/7/31",
             "2020/8/31","2020/9/30","2020/10/31","2020/11/30","2020/12/31",
             "2021/1/31","2021/2/28","2021/3/31","2021/4/30","2021/5/31","2021/6/30",
             "2021/7/31","2021/8/31","2021/9/30","2021/10/31","2021/11/30","2021/12/31"]

figname = ["2019-1","2019-2","2019-3","2019-4","2019-5","2019-6",
           "2019-7","2019-8","2019-9","2019-10","2019-11","2019-12",
          "2020-1","2020-2","2020-3","2020-4","2020-5","2020-6","2020-7",
           "2020-8","2020-9","2020-10","2020-11","2020-12","2021-1","2021-2","2021-3","2021-4","2021-5",
           "2021-6","2021-7","2021-8","2021-9","2021-10","2021-11","2021-12"]

for i in need_time[:28]:
    
    j = int(j)+1
    name = str(j)+"month_mean_"
  
##2020年5月之后只有12个站的数据
    if j>16:
        pos = pos2
        lons = pos2["经度"]
        lats = pos2["纬度"]
    else:
        pos = pos1
        lons = pos1["经度"]
        lats = pos1["纬度"]
    
    data  = pd.DataFrame(df[df["date"]==i]["O3"])
    min_data = data.min().values
    max_data = data.max().values
    
    #print(data)
    grid_lon = np.linspace(112.4, 113.8,140)
    grid_lat = np.linspace(34.7, 35.6,90)
    
    OK = OrdinaryKriging(lons, lats, data, variogram_model='linear',nlags=6)
    z1, ss1 = OK.execute('grid', grid_lon, grid_lat)
    #z1.shape
    
    xgrid, ygrid = np.meshgrid(grid_lon, grid_lat)
    df_grid = pd.DataFrame(dict(long=xgrid.flatten(),lat=ygrid.flatten()))
    df_grid["Krig_gaussian"] = z1.flatten()
    #df_grid.to_csv(name+".csv",encoding='utf_8_sig')
    
    sh = gpd.read_file(r'D:\arcgis\jiaozuo\final\jiaozuo_result.shp')
    #sh
    
    fig = plt.figure(figsize=(12, 8))
    ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
    ax.set_extent([112.5, 113.7, 34.75, 35.55], crs=ccrs.PlateCarree())
    
    province = shpreader.Reader(r"D:\arcgis\jiaozuo\final\jiaozuo_result.shp")
    ax.add_geometries(province.geometries(), crs=ccrs.PlateCarree(), linewidths=1,edgecolor='k',facecolor='none')
    #cdict=['white', '#A6F28F','#3DBA3D','#61B8FF','#0000FF','#FA00FA','#800040','#60281E']#自定义颜色列表
    cf = ax.contourf(xgrid, ygrid, z1, levels=np.linspace(120,130,21),cmap="Spectral_r", 
                     transform=ccrs.PlateCarree())
    cb = plt.colorbar(cf,shrink=0.8,pad=0.01)
    cb.set_label('O${_3}$($\mu$g/m${^3}$)',fontsize=15)
    
    shp2clip(cf, ax, r'D:\arcgis\jiaozuo\final\jiaozuo_result.shp')
    
    #ax.set_extent([112., 114., 34, 36])
    ax.scatter(lons,lats,c=data["O3"],s=90,ec="k",lw=0.5,cmap="Spectral_r",vmin = 120,vmax =130)

    ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=True))
    ax.yaxis.set_major_formatter(LatitudeFormatter())
    ax.set_xticks(np.arange(112.6,113.7,0.2))
    ax.set_yticks(np.arange(34.8,35.5,0.1))
    ax.tick_params(labelsize=18) #坐标轴刻度字体大小18
    
    ax.set_title(figname[j],fontsize=22)

    plt.show()
    fig.align_ylabels(ax)
    
    fig.savefig(name+'.png',bbox_inches = 'tight',dpi=300) 

为了严谨,coordinates_type="geographic"应该补充在OrdinaryKriging这个函数里。不过对于中纬度而言,加与不加差距很小了。克里金插值的coordinates_type参数 - Heywhale.com

相关文章

网友评论

      本文标题:如何根据站点污染物数据插值绘制出污染物的空间分布?⑨

      本文链接:https://www.haomeiwen.com/subject/slrnjrtx.html