Hi all,
I’ve been coding for a little over two years now, but it has mostly been to make plots and perform simple calculations. Most other times, I have colleagues who also code in Julia and have codes/functions/packages that simplified my coding experience. So I’m not that familiar with the technical jargon and I apologize beforehand if my explanation about my issue doesn’t make sense.
I’m trying to compile a database of hurricane-centered, passive microwave-based products from polar-orbiting satellites. Below is an example of one case, but I’m generalizing my code to apply to other cases. What I’ve done so far is converted the pixel location, originally in the form of latitude and longitude of the satellite pixel, to distance from the storm center. Below is an example of the surface precipitation field from Hurricane Matthew (2016) before and after I interpolated the distance.
Before I explain my issue, I have to first explain the data structure. The variables from the satellite product are given in terms of m x n dimension, where m corresponds to the cross track pixels and n corresponds to the along track scan. The cross-track pixels are perpendicular to the satellite trajectory as it traverses the earth, where the satellite scans along an arc. The along track scan is essentially how many cross-track scans are made, and is parallel to the trajectory of the satellite. All the variables from this dataset is given in this dimension, including the lat/lon (which I have converted to x and y having the same dimension).
So, if I want to plot the same data as above using PyPlot, but focusing on one or the other dimension, here’s what the code would look like:
# Remove missing/fill values
var[findin(var, -9999.0)] = NaN;
plt.subplot(1,2,1)
# Focus on cross-track pixels from index 20 to 22 for all along-track scan
plt.contourf(x[20:22, :], y[20:22, :], var[20:22, :], levels=(0:70))
plt.subplot(1,2,2)
# Focus on along-track scans from index 170:190 for all cross-track pixels
plt.contourf(x[:, 170:190], y[:, 170:190], var[:, 170:190], levels=(0:70))
Producing the plot:
The next step I would like to take is to interpolate the x and y values (again, distance relative to the storm center) onto a Cartesian grid from -800 km to 800 km along both x and y, with a constant grid spacing. I’ve looked at both the Julia Interpolations.jl package, and the Dierckx.jl package, but neither seem to be able to do what I want to do.
For the Interpolations package, to create an interpolation object, the input x and y values must each be a 1-d array, while the var values must have dimension [length(x), length(y)]. As you can see from the data structure that I am working with, that is not possible, since each index m and n in either x, y, or var is unique, because of the way the satellite scans.
I tried to use the Dierckx package, where a 2-d interpolation could be conducted if all x, y, and var have the same dimensions. So, I re-shaped my x, y, and var into 1-d arrays with length (m times n) and tried the Dierckx interpolation (copy and pasted from the examples given on the GitHub page). Below is what I have:
x_new = reshape(x, length(x))
y_new = reshape(y, length(y))
var_new = reshape(var, length(var))
spl = Spline2D(x_new, y_new, var_new; kx=3, ky=3, s=0.0)
But that crashed my Jupyter Notebook kernel (at least I think it did, 5 minutes went by without any indication that it finished running). I tried running the above code on a smaller subset of x, y, and var, and I got the error message, “The required storage space exceeds the available storage space: nxest or nyest too small, or s too small. Try increasing s.” and I’m not sure what that means. I increased s to 4000, but I still got the same answer. I have to admit that attempt was a desperate shot in the dark and I don’t really understand how the smoothing factor works.
On top of all this, the satellite data has missing values (-9999.0), and I’m not sure how the different interpolations package would handle that or NaNs. After interpolating the data onto a Cartesian grid with constant grid spacing, I would also like to fill the regions not covered by the satellite scan with fill values, so that I would have a 1600 x 1600 km grid with either a missing value or the actual value at each grid point.
Would anyone have an idea as to how I could go about doing this? I would greatly appreciate any help I can get.
