Untitled

avatar
unknown
plain_text
9 months ago
11 kB
6
Indexable
# 导入所需库
import time
import cartopy.crs as ccrs
# import cmaps
from osgeo import gdal
import matplotlib as mpl
import matplotlib.patches as mpatches
import matplotlib.pyplot as plt
import numpy as np
from cartopy.io.shapereader import Reader
from cartopy.mpl.ticker import LatitudeFormatter, LongitudeFormatter
from matplotlib.ticker import MultipleLocator

# from matplotlib.patheffects import withStroke

plt.rcParams.update({'font.size': 12})
# 添加指北针的函数
def add_north(ax, labelsize=10.5, loc_x=0.06, loc_y=0.94, width=0.04, height=0.09, pad=0.14):
    """
    画一个比例尺带'N'文字注释
    主要参数如下
    :param ax: 要画的坐标区域 Axes实例 plt.gca()获取即可
    :param labelsize: 显示'N'文字的大小
    :param loc_x: 以文字下部为中心的占整个ax横向比例
    :param loc_y: 以文字下部为中心的占整个ax纵向比例
    :param width: 指南针占ax比例宽度
    :param height: 指南针占ax比例高度
    :param pad: 文字符号占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='white')  # 将颜色更改为白色
    triangle = mpatches.Polygon([left, top, right, center], color='white', edgecolor='black')  # 添加黑色边框

    ax.text(s='N',
            x=minx + xlen*loc_x,
            y=miny + ylen*(loc_y - pad + height),
            fontsize=labelsize,
            horizontalalignment='center',
            verticalalignment='bottom',
            color='white')  # 文字颜色也更改为白色
    ax.add_patch(triangle)


# 添加比例尺的函数，单位是千米
'''
def add_scalebar(ax,lon0,lat0,length,h1 = 0.25,h2 = 0.15,h3 = 0.6,h4 = 0.05,size = 10.5):
    # h1，h2 主副刻度高度 h3刻度标签高度 h4 单位与比例尺的距离
    # style 3
    ax.hlines(y=lat0,  xmin = lon0, xmax = lon0+length/111, colors="black", ls="-", lw=1, label='%d km' % (length),zorder = 10)
    ax.vlines(x = lon0, ymin = lat0-0, ymax = lat0+h1, colors="black", ls="-", lw=1)
    ax.vlines(x = lon0+length/2/111, ymin = lat0-0, ymax = lat0+h2, colors="black", ls="-", lw=1)
    ax.vlines(x = lon0+length/111, ymin = lat0-0, ymax = lat0+h1, colors="black", ls="-", lw=1)
    ax.text(lon0+length/111,lat0+h3,'%d' % (length),horizontalalignment = 'center',fontsize = size)
    if length/2>=1:
        ax.text(lon0+length/2/111,lat0+h3,'%d' % (length/2),horizontalalignment = 'center',fontsize = size)
    else:
        ax.text(lon0+length/2/111,lat0+h3,'%.1f' % (length/2),horizontalalignment = 'center',fontsize = size)
    ax.text(lon0,lat0+h3,'0',horizontalalignment = 'center',fontsize = size)
    ax.text(lon0+length/111+h4,lat0+h3,'km',horizontalalignment = 'center',fontsize = size)
'''
# 添加比例尺的函数，单位是米

def add_scalebar(ax, lon0, lat0, length, h1=0.25, h2=0.15, h3=0.6, h4=0.05, size=10.5):
    # h1，h2 主副刻度高度 h3刻度标签高度 h4 单位与比例尺的距离
    length_m = length * 1000 # 将长度从公里转换为米
    ax.hlines(y=lat0, xmin=lon0, xmax=lon0+length_m/111000, colors="black", ls="-", lw=1, zorder=10) # 注意此处的111改为111000，因为地球上每度的距离约为111km或111000m
    ax.vlines(x=lon0, ymin=lat0, ymax=lat0+h1, colors="black", ls="-", lw=1)
    ax.vlines(x=lon0+length_m/2/111000, ymin=lat0, ymax=lat0+h2, colors="black", ls="-", lw=1) # 此处的length也要改
    ax.vlines(x=lon0+length_m/111000, ymin=lat0, ymax=lat0+h1, colors="black", ls="-", lw=1) # 此处的length也要改

    # 文本标签不再需要检查是否小于1，因为我们已经转换为米
    ax.text(lon0+length_m/111000, lat0+h3, '%d' % (length_m), horizontalalignment='center', fontsize=size)
    ax.text(lon0+length_m/2/111000, lat0+h3, '%d' % (length_m/2), horizontalalignment='center', fontsize=size)

    ax.text(lon0, lat0+h3, '0', horizontalalignment='center', fontsize=size)
    ax.text(lon0+length_m/111000+h4, lat0+h3, 'm', horizontalalignment='center', fontsize=size) # 单位改为'm'

def add_scalebar_w(ax,lon0,lat0,length,h1 = 0.25,h2 = 0.15,h3 = 0.6,h4 = 0.05,size = 10.5):
    # h1，h2 主副刻度高度,越大越长, h3刻度标签高度，越大数字越靠上 h4 单位与比例尺的距离, 越大越靠右
    # style 3
    ax.hlines(y=lat0,  xmin = lon0, xmax = lon0+length/111, colors="white", ls="-", lw=1, label='%d km' % (length),zorder = 10)
    ax.vlines(x = lon0, ymin = lat0-0, ymax = lat0+h1, colors="white", ls="-", lw=1)
    ax.vlines(x = lon0+length/2/111, ymin = lat0-0, ymax = lat0+h2, colors="white", ls="-", lw=1)
    ax.vlines(x = lon0+length/111, ymin = lat0-0, ymax = lat0+h1, colors="white", ls="-", lw=1)
    ax.text(lon0+length/111,lat0+h3,'%d' % (length*1000),horizontalalignment = 'center',fontsize = size, color="white")
    '''
    if length/2>=1:
        ax.text(lon0+length/2/111,lat0+h3,'%d' % (length*1000/2),horizontalalignment = 'center',fontsize = size, color="white")
    else:
        ax.text(lon0+length/2/111,lat0+h3,'%.1f' % (length*1000/2),horizontalalignment = 'center',fontsize = size, color="white")
    '''
    ax.text(lon0+length/2/111,lat0+h3,'%d' % (length*1000/2),horizontalalignment = 'center',fontsize = size, color="white")
    ax.text(lon0,lat0+h3,'0',horizontalalignment = 'center',fontsize = size, color="white")
    ax.text(lon0+length/111+h4,lat0+h3,'m',horizontalalignment = 'center',fontsize = size, color="white")



# outline = r'/mnt/c/Users/xue98/OneDrive - The University of Tokyo/研究生/2023-5-31 QGIS/dem_differential/Differential.tiff'
outline = r'dam.tif'
# outline = r'/mnt/c/Users/xue98/OneDrive - The University of Tokyo/研究生/2023-5-31 QGIS/Zonation of the Atami debris flow/motion_city_background.tiff'

proj=ccrs.PlateCarree()
fig = plt.figure(figsize=(8, 4.8), dpi = 200)

# region1 = [139.074608602, 139.0884550000000104,35.1097079999999977,35.1224270000000018]

# 定义绘图区域
# ax1 = fig.add_axes([0.05, 0.05 , 0.9, 0.9], projection=proj)
ax1 = fig.add_axes([0.04, 0.08 , 0.8, 0.84], projection=proj)
# ax1.xaxis.set_major_formatter(LongitudeFormatter(number_format='.3f'))
# ax1.yaxis.set_major_formatter(LatitudeFormatter(number_format='.3f'))
 
# 卫星影像
ds2 = gdal.Open(outline)
geotransform = ds2.GetGeoTransform()
xmin = geotransform[0]
ymax = geotransform[3]
xres = geotransform[1]
yres = geotransform[-1]
cols = ds2.RasterXSize
rows = ds2.RasterYSize
arr = ds2.ReadAsArray()
ds = None
xmax = xmin+xres*cols
ymin = ymax+yres*rows
extent=[xmin,xmax,ymin,ymax]

# arr = np.where(arr==0,255,arr)

arr = arr.transpose((1, 2, 0))
ax1.imshow(arr, extent=extent)
# 设置了extent右边就会留有空白
# ax1.set_extent(region1)


# 设置图例
legend_elements = [
    # mpatches.Patch(color='#2c7fb8', label='Source area', fill=True),
    plt.Line2D([], [], color='white', label='Cross Section', linestyle='--', linewidth=1.5),


    # plt.Line2D([], [], color='#b2182b', label='Aochi River', linewidth=1.5),
    # plt.Line2D([], [], color='#2166ac', label='Terayama River', linewidth=1.5),
    # plt.Line2D([], [], color='#67a9cf', label='Narusawa River', linewidth=1.5),
    
    # mpatches.Patch(color='#7fcdbb', label='Runout zone', fill=True),
    # mpatches.Patch(color='#edf8b1', label='Floodplain', fill=True),
    # mpatches.Patch(color='#fc8d62', label='Damaged construction', fill=True),
    # mpatches.Patch(color='#8da0cb', label='Undisturbed construction', fill=True),
    # mpatches.Patch(color='white', label='Railway', fill=False),

    
]
# legend = ax1.legend(handles=legend_elements, loc='lower left', ncol=2)
legend = ax1.legend(handles=legend_elements, loc='lower left')


# 设置刻度
ax1.xaxis.set_major_formatter(LatitudeFormatter(direction_label=True, degree_symbol='°', number_format='g', transform_precision=1e-08, dms=True, minute_symbol='′', second_symbol='″', seconds_number_format='.0f', auto_hide=True, decimal_point=None, cardinal_labels=None))
ax1.yaxis.set_major_formatter(LongitudeFormatter(direction_label=True, degree_symbol='°', number_format='g', transform_precision=1e-08, dms=True, minute_symbol='′', second_symbol='″', seconds_number_format='.0f', auto_hide=True, decimal_point=None, cardinal_labels=None))
# ax1.set_xticks(np.arange(region1[0] , region1[1]  + .001, .01))
# ax1.set_yticks(np.arange(region1[-2], region1[-1] + .001, .01))

# ax1.set_xticks([139.07384875000002, 139.0787175, 139.08358625])
# yticks = [35.11288775, 35.1160675, 35.11924725]
# ax1.set_yticks(yticks)

ax1.set_xticks([139.0751948728, 139.0757811435, 139.0763674143])
# yticks = [35.11924993, 35.119640777, 35.120031624]
yticks = [35.1193802123, 35.1199013417]
ax1.set_yticks(yticks)

y_major_formatter = ax1.yaxis.get_major_formatter()

ax1.set_yticklabels(yticks, rotation='vertical', verticalalignment='center')

ax1.yaxis.set_major_formatter(y_major_formatter)

# 在ax1上添加指北针和比例尺
# add_north(ax1)
add_scalebar_w(ax1, 139.075882, 35.120225, 0.1, h1=0.0025/35, h2=0.0015/35, h3=0.6/7000, h4=0.00012)
    
# 设置绘图区域的范围
# ax1.set_extent(region1, crs=proj)

# 在ax1上添加coastlines
ax1.coastlines()

# 添加注释annotate
plt.annotate('sediment control dam', 
             xy=(139.076069, 35.119663),  # 这个是数据坐标系中的点
             xytext=(0.45, 0.7),  # 这是轴域百分比坐标系中的点
             textcoords='figure fraction',  # 指定xytext的坐标系
             arrowprops=dict(facecolor='#74C476', 
                             alpha=0.8, 
                             arrowstyle='-|>', 
                             connectionstyle='arc3, rad=0.2',
                             lw=1),
             bbox=dict(boxstyle='round,pad=0.5', fc='white', ec='black', lw=1, alpha=0.8),
             color='black', 
             fontsize=12
            )
img = ax1.imshow(arr, extent=extent)

# Create a new axes for the colorbar
cax = fig.add_axes([0.83, 0.095, 0.03, 0.81])  # Adjust these values to move and resize the colorbar

# Create a colorbar
cmap = plt.get_cmap('RdBu')
# cmap = plt.get_cmap('RdBu_r')
norm = mpl.colors.Normalize(vmin=-10, vmax=10)
cb1 = mpl.colorbar.ColorbarBase(cax, cmap=cmap, norm=norm, orientation='vertical')
cb1.set_label('Elevation (m)')

# cbar = plt.colorbar(img, ax=ax1)
# cbar.set_label('Elevation (m)')
# cbar.set_ticks(np.linspace(-15, 15, num=6))


'''
# effects=[withStroke(linewidth=2.5, foreground="#0066ff")]
ax1.text(139.078813, 35.116576,'Check Dam', 
        color='red',  # 设置文字颜色
        fontsize=10,  # 设置文字大小
        # path_effects=effects,
        bbox=dict(facecolor='white', edgecolor='black'))
'''
# 保存图像
# plt.savefig('map.png', dpi=300, bbox_inches='tight')
plt.savefig('dam_label.png', dpi=300)
plt.show()
Editor is loading...
Leave a Comment