interpolated map using points python

Question

I need to obtain a colour map using interpolation of altitude records. I have a set of 34 points located along a region and in each point I have altitude record. I tried the code below, but it results in a map that cannot extrapolate to the State borders, besides the colourful map is not presentable (Fig1). I tried the np.meshgrid and did not work because I do not have the full area fillen in with points (it is not a perfect grid). I use 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

Is that possible to do with Python 3.7? Anyone could help me, please?

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.

I tried the plt.pcolormesh(lon,lat,shade,cmap=cmap,norm=norm,zorder=0), but didn't work since this function works for a perfect grid. Do you have any other suggestion? — bruce_skywalker, Aug 04 '20 at 13:02
I think @MadPhysicist's point is that what you are asking for is not interpolation, it looks like you have that covered. You are asking for extrapolation, which usually require to specify a few more hypothesis. — YvesQuemener, Jan 21 '22 at 11:36

pandas_as_pd · Answer 1 · 2022-02-21T22:49:27.987

To create an interpolation map, I recommend using cloupy

Firstly, the data should be cleaned up and properly structured to meet cloupy's data structure requirements for creating an interpolation map:

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')

The required data structure to create an interpolation map is the pandas.DataFrame object with the following column order: value, longitude, latitude. When the data meets the required data structure, we can create an interpolation map object:

import cloupy as cl
import numpy as np

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

The shapefile_path argument shows the path to the shapefile from which the state boundaries will be drawn. In this case, the shapefile from this website was used, but it could be any other state's shapefile (the state of Paraná in Brazil). The crs argument specifies the shapefile's coordinate system (usually epsg:4326, but it can be different). When the interpolation map object is properly prepared, we can specify drawing arguments and draw the map:

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),
)

This map is like the map in the question. The interpolation result is not satisfactory and the map style is poor. However, the imap.draw() method has many arguments which can be specified and change significantly the interpolation result or/and the map style. Exemplary map based on the same data but with different arguments set in the imap.draw() method:

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'
)

For more information on creating interpolation maps see cloupy's README or check cloupy's docstrings.

Note that

when cloupy extrapolates values, it assumes that the values outside the mesh grid do not decrease/increase - it tries to keep the values outside the mesh grid as the extreme values of the mesh grid, which should be a correct approach in case of altitude
the cloupy's map object (cloupy.m_MapInterpolation) does not always require passing a shapefile manually. The cloupy's map object allows you to draw boundaries from the default shapefile using the country argument. The default shapefile contains administrative boundaries from the entire world
check the add_shape argument in the draw() method. It allows you to combine the main shape (boundary shape) with additional shapes
the map content will always be cropped to the main shape (the boundary shape specified in the cloupy.m_MapInterpolation(...) object)

interpolated map using points python

1 Answers1

To create an interpolation map, I recommend using cloupy

Note that