2

I have measurements (PPI arc scans) taken with a doppler wind lidar. The data is stored in a pandas dataframe where rows represent azimuth angle and columns represent radial distance (input shape = 30x197). Link to example scan, (csv). I want to transform this to a cartesian coordinate system, and output a 2d array which is re-gridded into x,y coordinates instead of polar with the values stored in the appropriate grid cell. Interpolation (nearest neighbor) is ok and so is zero or NaN padding of areas where no data exists.

Ideally the X and Y grid should correspond to the actual distances between points, but right now I'm just trying to get this working. This shouldn’t be terribly difficult, but I’m having trouble obtaining the result I want. So far, I have working code which plots on a polar axis beautifully (example image) but this won't work for the next steps of my analysis.

I have tried many different approaches with scipy.interpolate.griddata, scipy.ndimage.geometric_transform, and scipy.ndimage.map_coordinates but haven't gotten the correct output. Here is an example of my recent attempt (df_polar is the csv file linked): 

# Generate polar and cartesian meshgrids
r = df_polar.columns
theta = df_polar.index
theta = np.deg2rad(theta)

# Polar meshgrid
rad_c, theta_c = np.meshgrid(r,theta)

# Cartesian meshgrid
X = rad_c * np.cos(theta_c)
Y = rad_c * np.sin(theta_c)
x,y = np.meshgrid(X,Y)

# Interpolate from polar to cartesian grid
new_grid = scipy.interpolate.griddata(
    (rad_c.flatten(), theta_c.flatten()), 
    np.array(df_polar).flatten(), (x,y), method='nearest')

The result is not correct at all, and from reading the documentation and examples I don't understand why. I would greatly appreciate any tips on where I have gone wrong. Thanks a lot!!

Elliot
  • 93
  • 1
  • 6
  • What shape is your `new_grid` at the end? I get `(5910, 5910)`. – Bill Apr 16 '18 at 05:08
  • I just realised a very similar question has been [answered before](https://stackoverflow.com/questions/47146855/translating-radial-data-to-a-cartesian-grid-for-surface-plot?rq=1). – Bill Apr 17 '18 at 01:25
  • Thanks for the link, despite a lot of searching I hadn't seen that post! – Elliot Apr 18 '18 at 17:11

1 Answers1

1

I think you might be feeding griddata the wrong points. It wants cartesian points and if you want the values interpolated over a regular x/y grid you need to create one and provide that too.

Try this and let me know if it produces the expected result. It's hard for me to tell if this is what it should produce:

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

df_polar = pd.read_csv('onescan.txt', index_col=0)

# Generate polar and cartesian meshgrids
r = pd.to_numeric(df_polar.columns)
theta = np.deg2rad(df_polar.index)

# Polar meshgrid
rad_c, theta_c = np.meshgrid(r, theta)

# Cartesian equivalents of polar co-ordinates
X = rad_c*np.cos(theta_c)
Y = rad_c*np.sin(theta_c)

# Cartesian (x/y) meshgrid
grid_spacing = 100.0   # You can change this
nx = (X.max() - X.min())/grid_spacing
ny = (Y.max() - Y.min())/grid_spacing
x = np.arange(X.min(), X.max() + grid_spacing, grid_spacing)
y = np.arange(Y.min(), Y.max() + grid_spacing, grid_spacing)
grid_x, grid_y = np.meshgrid(x, y)

# Interpolate from polar to cartesian grid
new_grid = griddata(
    (X.flatten(), Y.flatten()),
    df_polar.values.flatten(),
    (grid_x, grid_y),
    method='nearest'
)

The resulting values look something like this (with grid_spacing = 10 and flipping x and y):

import matplotlib.pyplot as plt
plt.imshow(new_grid.T, cmap='hot')

enter image description here

Clearly interpolate "nearest" needs taming...

Bill
  • 10,323
  • 10
  • 62
  • 85
  • Hi Bill, thanks so much for your help on this! I understand now that the (x,y) Cartesian meshgrid I was passing before was just a transformation of the polar points and needed to be re-gridded on a normally spaced plane. Thank you for discovering and correcting that. You are right that the nearest neighbor interpolation is acting strangely. If I re-run with linear interpolation (and also reverse the values in order to get the orientation right) then this is my result: https://i.imgur.com/LEh0xdS.png This looks great, but I should use the nearest neighbor method. Any tips going forward? – Elliot Apr 16 '18 at 18:25
  • Wow, that looks great. I can't offer more advice as I haven't used this before. Does the nearest neighbours offer any options such as limiting the distance of the NN? – Bill Apr 17 '18 at 01:12
  • If you really can't get `griddata` to do what you want you could do the nearest-neighbours interpolation yourself. [scipy.spatial.KDTree](https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.spatial.KDTree.html) or [sklearn.neighbors](http://scikit-learn.org/stable/modules/neighbors.html) might be useful. – Bill Apr 17 '18 at 01:21
  • Feel free to [acknowledge my answer](https://stackoverflow.com/help/someone-answers) @Elliot. – Bill Apr 17 '18 at 22:38
  • Thanks Bill, I have thought more about this and you are correct, I will need to do the NN interpolation manually as of the current state of scipy's griddata function. I will try to proceed by checking if the points are within the convex hull of the scan: https://stackoverflow.com/questions/16750618/whats-an-efficient-way-to-find-if-a-point-lies-in-the-convex-hull-of-a-point-cl – Elliot Apr 18 '18 at 17:15