0

I am trying to plot 2d terrain map with x,y and z (elevation). I followed the steps from the following link but I am getting very weird plot.

Python : 2d contour plot from 3 lists : x, y and rho?

I spent almost half day searching but got nowhere. 

import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate

# import data:
import xlrd
loc = "~/Desktop/Book4.xlsx"
wb = xlrd.open_workbook(loc)
sheet = wb.sheet_by_index(0)
sample=500

# Generate array:
x=np.array(sheet.col_values(0))[0:sample]
y=np.array(sheet.col_values(1))[0:sample]
z=np.hamming(sample)[0:sample][:,None]

# Set up a regular grid of interpolation points

xi, yi = np.meshgrid(x, y)

# Interpolate
rbf = scipy.interpolate.Rbf(x, y, z, function='cubic')
zi = rbf(xi, yi)
# Plot
plt.imshow(zi, vmin=z.min(), vmax=z.max(), origin='lower',
           extent=[x.min(), x.max(), y.min(), y.max()])
plt.colorbar()
plt.show()

The first of the following fig is what I am getting and the last one is how it should look like. This is what I am getting This is how it should look like

Any help shall be appreciated

Link to data file

gfdsal
  • 651
  • 1
  • 8
  • 25
  • if you have x,y,z columns you could use plt.tricontourf? – Derek Eden Oct 01 '19 at 23:29
  • Your data file has only 2 columns, _x_ and _y_ positions. But you need some values to interpolate... sounds like you are expecting elevation? Is there more data somewhere? (Related: I don't understand what the Hamming window comes into it for.) – Matt Hall Oct 02 '19 at 00:32
  • @kwinkunks I will have the elevation from another source, for now I am just simulating using Hamming window the elevation values. Regardless in both the cases my elevation input is a 1D vector of same length as the other two. – gfdsal Oct 02 '19 at 11:24
  • @DerekEden plt.tricontourf gives weird triangulated plot which I am not targetting. I was the diffused plot as generated in second picture. If I do plain scatter plot of x,y I do get the shape of the second plot but without the density colorings. – gfdsal Oct 02 '19 at 11:27
  • I think the problem is that you're not really simulating anything meaningful because of the order of the points. So the result is not at all smooth and the interpolation doesn't work (it ends up with a singular matrix). You can relax the condition that the interpolation goes through all the data points by using the `smooth` argument in the interpolator. I'll post an example. – Matt Hall Oct 02 '19 at 13:00

1 Answers1

1

I think the problem is that the data you're giving it is not smooth enough to interpolate with the default parameters. Here's one approach, using mgrid instead of meshgrid:

import numpy as np
import pandas as pd
from scipy.interpolate import Rbf

# fname is your data, but as a CSV file.
data = pd.read_csv(fname).values
x, y = data.T

x_min, x_max = np.amin(x), np.amax(x)
y_min, y_max = np.amin(y), np.amax(y)

# Make a grid with spacing 0.002.
grid_x, grid_y = np.mgrid[x_min:x_max:0.002, y_min:y_max:0.002]

# Make up a Z.
z = np.hamming(x.size)

# Make an n-dimensional interpolator.
rbfi = Rbf(x, y, z, smooth=2)

# Predict on the regular grid.
di = rbfi(grid_x, grid_y)

Then you can look at the result:

import matplotlib.pyplot as plt

plt.imshow(di)

I get:

the result of the interpolation

I wrote a Jupyter Notebook on this topic recently, check it out for a few other interpolation methods, like kriging and spline fitting.

Matt Hall
  • 7,614
  • 1
  • 23
  • 36
  • I will review the notebook. I tried to reproduce your results. I am getting the same image but quite small in size and pixelated. What parameters can I adjust to obtain the proper resolution as above. Secondly the axes in the reproduced plot are absolute including origin and not relative as in the original data. – gfdsal Oct 02 '19 at 14:34
  • 1
    You can increase the size of the plot by starting the plot with `plt.figure(figsize=(10, 15))` (for example). The pixellation is a function of the grid step size, which I gave as 0.002 -- but you can use whatever you like. Make it smaller for more grid cells. You can also use the `interpolation` argument in `imshow` to smoth the display (but not the data). – Matt Hall Oct 02 '19 at 14:37
  • This makes sense. Also I was looking for how can I obtain output Z by querying a specific (x,y) value that is not in the array of data (x,y,z). I believe your interpolation method might be helpful in that regard as well as I Bayesian optimization that will query for specific values based on acquisition function. So to run an experiment if the specific values are not in the data, I have to do some interpolation to get those values to be fed into BO which uses GP. – gfdsal Oct 02 '19 at 14:44
  • 1
    The `rbfi()` interpolation function can look up (x, y) locations. Other methods have similar functions, eg `predict()` on an `sklearn` model. Good luck! – Matt Hall Oct 02 '19 at 14:48