import matplotlib.pyplot as plt
import numpy as np
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import pandas as pd
import maskout
import matplotlib.patches as mpatches
from cartopy.io.shapereader import Reader
from matplotlib import rcParams
config = {"font.family":'Times New Roman',"font.size": 16,"mathtext.fontset":'stix'}
rcParams.update(config)
df1 = pd.read_excel(r"F:/Rpython/lp37/henanmap/henan.xlsx")
lat = df1['lat']
lon = df1['lon']
PM25 = df1['PM25']
# 設(shè)置經(jīng)緯度并用scipy包中的函數(shù)插值
from scipy.interpolate import Rbf
olon = np.linspace(110, 117, 100)
olat = np.linspace(31, 37, 100)
olon, olat = np.meshgrid(olon, olat)
# 插值處理,‘linear'，‘nearest'，‘cubic'
# cubic, gaussian, inverse_multiquadric, linear, multiquadric, quintic, thin_plate
# rain_data_new = griddata((lon,lat), data, (olon,olat), method='linear')
func1=Rbf(lon,lat,PM25,function='linear')
PM25=func1(olon,olat)
print(olon)
print(olat)
print(PM25)
print(olat.shape)
print(PM25.shape)
# 建立畫布
fig2 = plt.figure(figsize=(16, 12))
proj = ccrs.PlateCarree()
leftlon, rightlon, lowerlat, upperlat = (110, 117, 31, 37)  # 根據(jù)上下限確定范圍，至少為10°
# 春 ------------------------------------------------------------------------------------
#左 底 寬 高
ax=fig2.add_subplot(1,1,1,projection=ccrs.PlateCarree())
# 在畫布的絕對(duì)坐標(biāo)建立子圖
ax.set_extent([leftlon, rightlon, lowerlat, upperlat],crs=ccrs.PlateCarree())
# 設(shè)置地圖屬性:加載國(guó)界、海岸線，land 為灰色
ax.stock_img()
# 繪制河南省行政邊界
ticks=np.arange(30,80,5)
cf=ax.contourf(olon,olat,PM25,levels=ticks,cmap='gist_rainbow',transform=ccrs.PlateCarree(),extend='both')
clip1=maskout.shp2clip(cf,ax,'F:/Rpython/lp37/henanmap/henan1')
cf=plt.colorbar(cf,ticks=ticks,shrink=0.96,orientation='vertical',extend='both',pad=0.01,aspect=35)
ax.add_feature(cfeature.COASTLINE.with_scale('50m'),linewidth=0.5,zorder=2,color='k')# 添加海岸線
ax.add_feature(cfeature.LAKES.with_scale('50m'))
ax.add_feature(cfeature.RIVERS.with_scale('50m'))
ax.add_feature(cfeature.OCEAN.with_scale('50m'))
ax.add_feature(cfeature.LAND.with_scale('50m'))
ax.add_geometries(Reader(r'F:/Rpython/lp37/henanmap/henan1.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=1.2,zorder=1)
ax.add_geometries(Reader(r'F:/Rpython/lp27/data/river1.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='b',linewidth=0.2)
ax.add_geometries(Reader(r'F:/Rpython/lp27/data/china1.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='r',linewidth=1.2)
ax.add_geometries(Reader(r'F:/Rpython/lp27/data/china2.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='r',linewidth=0.8)
ax.add_geometries(Reader(r'F:/Rpython/lp27/data/ne_50m_lakes.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=0.2)
ax.add_geometries(Reader(r'F:/Rpython/lp37/data37/river/1級(jí)河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='RoyalBlue',linewidth=0.4)
ax.add_geometries(Reader(r'F:/Rpython/lp37/data37/river/2級(jí)河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='DodgerBlue',linewidth=0.3)
ax.add_geometries(Reader(r'F:/Rpython/lp37/data37/river/3級(jí)河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='DeepSkyBlue',linewidth=0.2)
ax.add_geometries(Reader(r'F:/Rpython/lp37/data37/river/4級(jí)河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='SkyBlue',linewidth=0.15)
ax.add_geometries(Reader(r'F:/Rpython/lp37/data37/river/5級(jí)河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='LightSkyBlue',linewidth=0.05)
ax.add_geometries(Reader(r'F:/Rpython/lp37/data37/river/主要湖泊.shp').geometries(),ccrs.PlateCarree(),edgecolor='none',linewidth=0,facecolor='#BEE8FF')
# 以下6條語(yǔ)句是定義地理坐標(biāo)標(biāo)簽格式
ax.set_xticks(np.arange(leftlon, rightlon+0.1, 1),crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(lowerlat, upperlat+0.1, 1),crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter()
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
plt.savefig('F:/Rpython/lp37/plot185.4.png',dpi=600,bbox_inches='tight',pad_inches=0)
plt.show()

# 添加比例尺，指南針
def add_north(ax, labelsize=18, loc_x=0.95, loc_y=0.99, width=0.06, height=0.09, pad=0.14):
    """
    畫一個(gè)比例尺帶'N'文字注釋
    主要參數(shù)如下
    :param ax: 要畫的坐標(biāo)區(qū)域 Axes實(shí)例 plt.gca()獲取即可
    :param labelsize: 顯示'N'文字的大小
    :param loc_x: 以文字下部為中心的占整個(gè)ax橫向比例
    :param loc_y: 以文字下部為中心的占整個(gè)ax縱向比例
    :param width: 指南針占ax比例寬度
    :param height: 指南針占ax比例高度
    :param pad: 文字符號(hào)占ax比例間隙
    :return: None
    """
    minx, maxx = ax.get_xlim()
    miny, maxy = ax.get_ylim()
    ylen = maxy - miny
    xlen = maxx - minx
    left = [minx + xlen*(loc_x - width*.5), miny + ylen*(loc_y - pad)]
    right = [minx + xlen*(loc_x + width*.5), miny + ylen*(loc_y - pad)]
    top = [minx + xlen*loc_x, miny + ylen*(loc_y - pad + height)]
    center = [minx + xlen*loc_x, left[1] + (top[1] - left[1])*.4]
    triangle = mpatches.Polygon([left, top, right, center], color='k')
    ax.text(s='N',
            x=minx + xlen*loc_x,
            y=miny + ylen*(loc_y - pad + height),
            fontsize=labelsize,
            horizontalalignment='center',
            verticalalignment='bottom')
    ax.add_patch(triangle)
#-----------函數(shù)：添加比例尺--------------
def add_scalebar(ax,lon0,lat0,length,size=0.01):
    '''
    ax: 坐標(biāo)軸
    lon0: 經(jīng)度
    lat0: 緯度
    length: 長(zhǎng)度
    size: 控制粗細(xì)和距離的
    '''
    # style 3
    ax.hlines(y=lat0,  xmin = lon0, xmax = lon0+length/111, colors="black", ls="-", lw=1,linewidth=3, label='%d km' % (length))
    ax.vlines(x = lon0, ymin = lat0-size, ymax = lat0+size, colors="black", ls="-", lw=1,linewidth=3)
    ax.vlines(x = lon0+length/2/111, ymin = lat0-size, ymax = lat0+size, colors="black", ls="-", lw=1,linewidth=3)
    ax.vlines(x = lon0+length/111, ymin = lat0-size, ymax = lat0+size, colors="black", ls="-", lw=1,linewidth=3)
    ax.text(lon0+length/111,lat0+size+0.05,'%d' % (length),horizontalalignment = 'center')
    ax.text(lon0+length/2/111,lat0+size+0.05,'%d' % (length/2),horizontalalignment = 'center')
    ax.text(lon0,lat0+size+0.05,'0',horizontalalignment = 'center')
    ax.text(112.28,31.45,'km',horizontalalignment = 'center')

def drow_the_scale(y,x,text,length = 1.5,lw = 5):
    #畫比例尺函數(shù)
    # y代表比例尺所在緯度
    # x代表比例尺開始的經(jīng)度
    # text代表比例尺最后刻度值
    # length代表比例尺的長(zhǎng)度，單位為多少個(gè)經(jīng)度
    # lw代表比例尺的寬度
    step = length/5#計(jì)算步長(zhǎng)，畫五格
    #畫黑白線五條
    plt.hlines(y=y,xmin=x,xmax=x + step,colors="black", ls="-", lw=lw)
    plt.hlines(y=y,xmin=x + step,xmax=x + step*2,colors="white", ls="-", lw=lw)
    plt.hlines(y=y,xmin=x + step*2,xmax=x + step*3,colors="black", ls="-", lw=lw)
    plt.hlines(y=y,xmin=x + step*3,xmax=x + step*4,colors="white", ls="-", lw=lw)
    plt.hlines(y=y,xmin=x + step*4,xmax=x + step*5,colors="black", ls="-", lw=lw)
    #畫長(zhǎng)刻度兩個(gè)
    plt.vlines(x = x, ymin = y - (lw/100) *3, ymax = y + lw/100, colors="black", ls="-", lw=1)
    plt.vlines(x = x + length, ymin = y - (lw/100) *3, ymax = y + lw/100, colors="black", ls="-", lw=1)
    #畫段刻度四個(gè)
    plt.vlines(x = x + step, ymin = y - (lw/100) *2, ymax = y + lw/100, colors="black", ls="-", lw=1)
    plt.vlines(x = x + step*2, ymin = y - (lw/100) *2, ymax = y + lw/100, colors="black", ls="-", lw=1)
    plt.vlines(x = x + step*3, ymin = y - (lw/100) *2, ymax = y + lw/100, colors="black", ls="-", lw=1)
    plt.vlines(x = x + step*4, ymin = y - (lw/100) *2, ymax = y + lw/100, colors="black", ls="-", lw=1)
    #寫字，0，500，km
    plt.text(x,y - (lw/100) *7,'0',horizontalalignment = 'center')
    plt.text(x + length,y - (lw/100) *7,text,horizontalalignment = 'center')
    plt.text(x + length/2,y + (lw/100)*2,'km',horizontalalignment = 'center')
def drowscale(extent,scale_y,scale_x,scale_text,step = 5,lw = 10,scale_length = 1.5,scale_lw = 5):
    # 畫地圖黑白邊框和比例尺
    # extent：表示四周經(jīng)緯度[west, east, south, north]
    # scale_y,scale_x,scale_text：代表比例尺的位置，緯度，經(jīng)度，刻度值
    # step：表示步長(zhǎng)，一格代表幾個(gè)經(jīng)緯度
    # lw：代表邊框?qū)挾?
    # scale_length：代表比例尺長(zhǎng)度（單位為經(jīng)度例如1.5個(gè)經(jīng)度）
    # scale_lw：代表比例尺寬度
    for y in [extent[2],extent[3]] :#畫上下兩邊框
        xmin = extent[0]
        while (xmin < extent[1]):
            plt.hlines(y=y,xmin=xmin,xmax=xmin+step,colors="white", ls="-", lw=lw)
            xmin = xmin+step*2
        xmin = extent[0]+step
        while (xmin < extent[1]):
            plt.hlines(y=y,xmin=xmin,xmax=xmin+step,colors="black", ls="-", lw=lw)
            xmin = xmin+step*2
    for x in [extent[0],extent[1]] :#畫左右兩邊狂
        ymin = extent[2]
        while (ymin < extent[3]):
            plt.vlines(x = x, ymin = ymin, ymax = ymin+step, colors="black", ls="-", lw=lw)
            ymin = ymin+step*2
        ymin = extent[2]+step
        while (ymin < extent[3]):
            plt.vlines(x = x, ymin = ymin, ymax = ymin+step, colors="white", ls="-", lw=lw)
            ymin = ymin+step*2
    drow_the_scale(scale_y,scale_x,scale_text)#畫比例尺