0

我需要使用高度记录的插值来获得彩色地图。我有一组位于一个区域的 34 个点，每个点都有海拔记录。我尝试了下面的代码，但它导致地图无法推断到州边界，此外彩色地图无法呈现（图 1）。我尝试了 np.meshgrid 并没有工作，因为我没有用点填充整个区域（它不是一个完美的网格）。我使用 Python 3.7。

在此处输入图像描述

#    Lat     Lon    Altitude
1   -25.53  -48.51  4.50
2   -25.30  -48.49  59.00
3   -25.16  -48.32  40.00
4   -25.43  -49.26  923.50
5   -24.49  -49.15  360.00
6   -25.47  -49.46  910.00
7   -23.30  -49.57  512.00
8   -25.13  -50.01  880.00
9   -23.00  -50.02  450.00
10  -23.06  -50.21  440.00
11  -24.78  -50.00  1008.80
12  -25.27  -50.35  893.00
13  -24.20  -50.37  768.00
14  -23.16  -51.01  484.00
15  -22.57  -51.12  600.00
16  -23.54  -51.13  1020.00
17  -25.21  -51.30  1058.00
18  -23.30  -51.32  746.00
19  -26.29  -51.59  1100.00
20  -26.25  -52.21  930.00
21  -25.25  -52.25  880.00
22  -23.31  -51.13  566.00
23  -23.40  -51.91  542.00
24  -23.05  -52.26  480.00
25  -24.40  -52.34  540.00
26  -24.05  -52.36  616.40
27  -23.40  -52.35  530.00
28  -26.07  -52.41  700.00
29  -25.31  -53.01  513.00
30  -26.05  -53.04  650.00
31  -23.44  -53.17  480.00
32  -24.53  -53.33  660.00
33  -25.42  -53.47  400.00
34  -24.18  -53.55  310.00

这可能与Python 3.7有关吗？任何人都可以帮助我，好吗？

import numpy as np
import matplotlib.pyplot as plt  
from matplotlib.colors import BoundaryNorm  #Organizes the color scale
import pandas as pd

import os
os.environ['PROJ_LIB'] = r'C:\Users\User\Anaconda3\pkgs\proj4-5.2.0-h6538335_1006\Library\share'

from mpl_toolkits.basemap import Basemap    #plots the maps in geographical coordinates

#Open a graphic window
plt.figure(figsize=(7.3, 4))
plt.subplots_adjust(left=.06,right=.99,top=.99,bottom=.05)

longlat=pd.read_excel(r'C:/Users/User/All data.xlsx')

lat=longlat.iloc[:,1].values
lon=longlat.iloc[:,2].values
div=1

#Prepare the lon and lat with the area to be ploted (acording to a lon and lat array)
mny=np.floor(min(lat)/div)*div; mxy=np.ceil(max(lat)/div)*div
sty=10 if (mxy-mny < 60) else 15
if(mxy-mny < 20): sty=4
if(mxy-mny > 100): sty=30

circles=np.arange(mny,mxy+sty,sty).tolist()
mnx=np.floor(min(lon)/div)*div; mxx=np.ceil(max(lon)/div)*div
stx=15 if (mxx-mnx < 100) else 30
if(mxx-mnx > 300): stx=45
if(mxx-mnx < 20): stx=4
meridians=np.arange(mnx,mxx+stx,stx).tolist()

mm = Basemap(llcrnrlon=mnx,llcrnrlat=mny,urcrnrlon=mxx,\
            urcrnrlat=mxy,resolution='l',area_thresh=10000.,\
            projection='cyl')     #This sets the coordinates in space, but won't plot any line
mm.drawcountries(zorder=1);mm.drawcoastlines(zorder=1,linewidth=.5)  #Continents will be ploted here
paral=mm.drawparallels(circles,linewidth='0.1',labels=[1,0,0,0],fontsize=8)  #Draw paralels

for m in paral: 
    try: 
        paral[m][1][0].set_rotation(90)
    except:
        dummy=""
mm.drawmeridians(meridians[1:],linewidth='0.1',labels=[0,0,0,1],fontsize=8)  #Draw meridians
plt.box(on=None)  #No frame around the map

inshp='C:/Users/User/Desktop/BR_States/estados_br_pol.shp'  #Input shp 

from shapefileBR_PR import shapefile_BR
states=shapefile_BR(inshp)

for i in states.keys(): plt.plot(states[i][0,:],states[i][1,:],'k',linewidth=0.5) #Plots states

lvl=(0,1200,100) #min, max and step of the shade
#Defining your color scheme
cmap=plt.get_cmap('RdBu')  #RdBu is a scale varing from red to blue. The opposite (blue to red) use 'RdBu_r'
levels=np.arange(lvl[0],lvl[1]+lvl[2],lvl[2])  #Define the levels (or range) you want to plot
norm=BoundaryNorm(levels,ncolors=cmap.N,clip=True)

shade = longlat.iloc[:,3].values

#Plotting values in irregular grid as filled countours
ct=plt.tricontourf(lon,lat,shade,levels=levels,cmap=cmap,norm=norm,zorder=0)  #shade should be an array [nlat,nlon]. Here lon and lat are arrays with the longitude and latitude of each station 
plt.scatter(lon,lat,marker='o',color='k',s=20,zorder=10) #this will plot the location of your stations as circles.
4

1 回答 1

0

要创建插值图，我建议使用cloupy

首先，数据应该被清理并适当地结构化以满足cloupy创建插值图的数据结构要求：

data = '''
1   -25.53  -48.51  4.50
2   -25.30  -48.49  59.00
3   -25.16  -48.32  40.00
4   -25.43  -49.26  923.50
5   -24.49  -49.15  360.00
6   -25.47  -49.46  910.00
7   -23.30  -49.57  512.00
8   -25.13  -50.01  880.00
9   -23.00  -50.02  450.00
10  -23.06  -50.21  440.00
11  -24.78  -50.00  1008.80
12  -25.27  -50.35  893.00
13  -24.20  -50.37  768.00
14  -23.16  -51.01  484.00
15  -22.57  -51.12  600.00
16  -23.54  -51.13  1020.00
17  -25.21  -51.30  1058.00
18  -23.30  -51.32  746.00
19  -26.29  -51.59  1100.00
20  -26.25  -52.21  930.00
21  -25.25  -52.25  880.00
22  -23.31  -51.13  566.00
23  -23.40  -51.91  542.00
24  -23.05  -52.26  480.00
25  -24.40  -52.34  540.00
26  -24.05  -52.36  616.40
27  -23.40  -52.35  530.00
28  -26.07  -52.41  700.00
29  -25.31  -53.01  513.00
30  -26.05  -53.04  650.00
31  -23.44  -53.17  480.00
32  -24.53  -53.33  660.00
33  -25.42  -53.47  400.00
34  -24.18  -53.55  310.00
'''

clean_data = []
for row in data.split('\n'):  
    row = row.split(' ')
    row = [value for value in row if value != '']
    clean_data.append(row)
clean_data = clean_data[1:-2] # delete empty values

import pandas as pd
df = pd.DataFrame(clean_data)
df = df.drop(0, axis=1)
df = df.iloc[:, [2, 1, 0]] # meet cloupy's data structure requirements

df = df.astype('float64')

创建插值图所需的数据结构是具有以下列顺序的 pandas.DataFrame 对象：值、经度、纬度。当数据满足所需的数据结构时，我们可以创建插值图对象：

import cloupy as cl
import numpy as np

imap = cl.m_MapInterpolation(
    shapefile_path='shapefiles/41UFE250GC_SIR.shp', 
    crs='epsg:4326', 
    dataframe=df
)

该shapefile_path参数显示将绘制状态边界的 shapefile 的路径。在这种情况下，使用了来自本网站的 shapefile ，但它可以是任何其他州的 shapefile（巴西的巴拉那州）。参数指定 shapefile的crs坐标系（通常是epsg:4326，但可以不同）。当插值地图对象准备好后，我们可以指定绘图参数并绘制地图：

imap.draw(
    levels=np.arange(0, 1200, 100),
    cmap='RdBu',
    interpolation_method='linear',
    interpolation_within_levels=True, # do not return values below 0, because the elevation must be equal to or greater than 0 
    boundaries_lw=0.2,
    show_points=True,
    show_grid=True,
    xticks=[-50],
    yticks=[-23, -27],
    figsize=(5, 4),
)

map_v1

这张地图就像问题中的地图。插值效果不理想，地图风格较差。但是，该imap.draw()方法有许多参数可以指定并显着改变插值结果或/和地图样式。基于相同数据但在imap.draw()方法中设置不同参数的示例地图：

imap.draw(
    levels=np.arange(0, 1550, 50),
    cmap='terrain',
    interpolation_method='cubic', # default value
    interpolation_within_levels=True,
    boundaries_lw=0.2,
    show_points=False, # default value
    show_grid=False, # default value
    xticks=np.arange(-55, -47, 1),
    yticks=np.arange(-27, -21, 1),
    figsize=(5, 4),
    show_contours=True,
    contours_levels=np.arange(250, 1501, 250),
    clabels_add=[
        (-53.6, -23.7), # specify manually the position of contour labels
        (-53.5, -25),
        (-52.5, -25.5),
        (-51, -25.5),
        (-51, -24.3),
    ],
    cbar_ticks=np.arange(0, 1550, 250),
    cbar_title='Elevation [m]',
    title='Topography in the state of Paraná (Brazil)',
    title_x_position=0.5,
    title_ha='center'
)

map_v2

有关创建插值图的更多信息，请参阅cloupy 的 README或检查 cloupy 的文档字符串。

注意

  • 当 cloupy 外推值时，它假定网格外的值不会减少/增加 - 它试图将网格外的值保持为网格的极值，这在海拔高度的情况下应该是正确的方法

  • cloupy 的地图对象 ( cloupy.m_MapInterpolation) 并不总是需要手动传递 shapefile。cloupy 的 map 对象允许您使用country参数从默认 shapefile 绘制边界。默认 shapefile 包含来自整个世界的管理边界

  • 检查方法add_shape中的参数draw()。它允许您将主要形状（边界形状）与其他形状结合起来

  • 地图内容将始终被裁剪为主要形状（cloupy.m_MapInterpolation(...)对象中指定的边界形状）

于 2022-02-21T03:38:02.957 回答