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.

Answer 1

Edit – Updated Answer

The cloupy package is not supported anymore due to outdated functionalities & unmaintainable code. See the README in the package repository for more information.

However, there is a new tool that perfectly addresses this problem – the GeoKrige package. This section of the documentation is especially helpful when it comes to creating interpolation maps in Python.

---

Creating interpolation map with GeoKrige

The GeoKrige package enables interpolation of data using Kriging Methods. To gain a better understanding of the GeoKrige package and Kriging Methods themselves, please refer to the GeoKrige documentation. A good starting point is this section.

1. Transform the data

Transform the data into numpy.ndarray objects, so that the GeoKrige package can operate on it:

import numpy as np

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

lines = data.strip().splitlines()
values = [[float(x) for x in line.split()[1:]] for line in lines]
data = np.array(values)

X = np.column_stack([data[:, 1], data[:, 0]])
y = data[:, 2]

2. Load the data & prepare the Kriging Model for interpolation

from geokrige.methods import OrdinaryKriging

kgn = OrdinaryKriging()
kgn.load(X, y)

kgn.variogram(plot=False)
kgn.fit(model='exp', plot=False)

Now, the Kriging Model is ready to predict values for the points with unknown values. The Kriging Model can be adjusted in many ways to meet the user's expectations, but the default settings in GeoKrige also handle many cases well. For more information about different components of the Kriging Model, please refer to the documentation.

3. Create a mesh grid on which predictions will be performed

import geopandas as gpd
from geokrige.tools import TransformerGDF

shp_file = 'shapefile/qj614bt0216.shp'
prediction_gdf = gpd.read_file(shp_file).to_crs(crs='EPSG:4326')

transformer = TransformerGDF()
transformer.load(prediction_gdf)

meshgrid = transformer.meshgrid(density=2)
mask = transformer.mask()

The shapefile has been downloaded from here. It has been read with the GeoPandas package and stored in a GeoDataFrame object.

The GeoKrige package facilitates an easy transformation of shapefiles into mesh grids on which interpolation is usually performed. Additionally, it is also possible to create effortlessly a masking array, which can be later handy in the visualization process.

Now the meshgrid object can be passed to the predict method of the Kriging Model.

4. Predict values & mask them to the shapefile boundaries

X, Y = meshgrid
Z = kgn.predict(meshgrid)

Z[~mask] = None

Here, the variables that will be used in the visualization are created. The mesh grid is divided into separate X and Y matrices, and another matrix with interpolated values is created – Z.

Note that the mask object is a boolean matrix – a True value signifies that a given point resides within the boundaries of the loaded shapefile. Thus, the masking matrix is inverted and the None values are assigned to the points that are located outside the boundaries. In this manner, the matplotlib package will make appropriate pixels empty.

5. Visualize the interpolation results

import matplotlib.pyplot as plt

fig, ax = plt.subplots()

prediction_gdf.plot(facecolor='none', edgecolor='black', linewidth=1.5, zorder=5, ax=ax)
cbar = ax.contourf(X, Y, Z, cmap='terrain', levels=np.arange(0, 1300, 50), extend='min')

# Cbar aligned with the plot on the right side
cax = fig.add_axes([0.93, 0.134, 0.02, 0.72])
colorbar = plt.colorbar(cbar, cax=cax, orientation='vertical')

clabels = ax.contour(X, Y, Z, levels=np.arange(250, 1501, 200), colors='k', linewidths=0.5)
ax.clabel(clabels, fontsize=8)

ax.grid(lw=0.2)
ax.set_title('Elevation of Parana State (Brazil)', fontweight='bold', pad=15)
ax.set_xlim(-55, -47.7)
ax.set_ylim(-27, -22.2)

plt.show()

---

Outdated Answer

Please, see the first two paragraphs of this answer – the cloupy package is no longer supported and there are better tools to resolve this problem.

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)