Search code examples
pythoninterpolationcartopycontourfdiscrete

Spatial interpolation of discrete points onto x/y coordinate mesh grid in Python


I'm still very new to programming and trying to create a contour plot of alkalinity across Hawaii using Cartopy. I will need to interpolate the point values called MODIFIED_TA against an x-y mesh grid but have not been able to figure out how to do this. The code I'm using is:

import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import cartopy.crs as ccrs
import cartopy.mpl.ticker as cticker
import statistics
from scipy.interpolate import UnivariateSpline
import numpy as np
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import warnings

warnings.filterwarnings("ignore") # ignoring the warning prompts.

%matplotlib inline

fig = plt.figure(figsize=(15,15))

ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree(central_longitude=-170))

landgreen = cfeature.NaturalEarthFeature('physical', 'land', '110m', 
                                        edgecolor='face', facecolor='green')


oceanaqua = cfeature.NaturalEarthFeature('physical', 'ocean', '110m', 
                                       edgecolor='face', facecolor='aqua')


ax.set_extent([-151.5, -162, 18, 24], ccrs.PlateCarree())
ax.set_title('TOTAL ALKALINITY')
ax.add_feature(landgreen)
ax.add_feature(cfeature.OCEAN, color = 'k')
ax.gridlines(draw_labels=True)


lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.grid(linewidth=2, color='black', alpha=0.5, linestyle='--')


lons = all_data.LONGITUDE[:]
lats = all_data.LATITUDE[:]
alk = all_data.MODIFIED_TA[:]
x1,y1 = np.meshgrid(lons,lats)
z1,z2 = np.meshgrid(all_data.MODIFIED_TA,all_data.MODIFIED_TA)


plt.tricontourf(lons,lats,alk, transform=ccrs.PlateCarree(), cmap=cm.gist_rainbow)
plt.colorbar(shrink=0.5)
plt.title('$A_{T}$ VALUES', color = 'k', fontname='Times New Roman',size = 23)

plt.plot()

The result is nothing like what I was hoping for and again, I'm not sure how to interpolate this value so that it comes out as a smooth gradient across the x/y coordinate grid. Any help would be greatly appreciated!

See output here


Solution

  • It's hard to tell for sure without being able to see your data. I tried to create a MRE and it worked. I would start by seeing if this works.

    import cartopy.feature as cfeature
    import matplotlib.pyplot as plt
    import cartopy.crs as ccrs
    import cartopy.mpl.ticker as cticker
    import numpy as np
    from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
    
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree(central_longitude=-170))
    
    ax.set_extent([-151.5, -162, 18, 24], ccrs.PlateCarree())
    ax.add_feature(cfeature.OCEAN)
    ax.gridlines(draw_labels=True)
    
    
    lon_formatter = cticker.LongitudeFormatter()
    lat_formatter = cticker.LatitudeFormatter()
    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)
    ax.grid(linewidth=2, color='black', alpha=0.5, linestyle='--')
    
    lons = np.random.random(80) * 7 - 160
    lats = np.random.random(80) * 4 + 19
    alk = np.cos(lons * 10 * np.pi / 180) * np.sin(lats * 20 / 180)
    
    plt.plot(lons, lats, 'k.', transform = ccrs.PlateCarree())
    plt.tricontourf(lons,lats,alk, transform=ccrs.PlateCarree(), alpha = 0.5)
    plt.colorbar(shrink=0.5)
    plt.title('$A_{T}$ VALUES', color = 'k', fontname='Times New Roman',size = 23)
    

    map

    If it does work, then what I'd look at would include:

    • What are the dimensions of all_data.LONGITUDE, all_data.LATITUDE, all_data.MOTIFIED_TA?
    • Are there duplicate values?
    • Does it work when you plot it outside of a projection?

    If my example does not work, then it suggests there is something about your install, in which case update it if you can. If the problem still persists, perhaps there is a bug cartopy that needs reporting or a conflict with other packages.

    Sorry, I cannot help further.