How to extract a One dimensional profile along a curve using Geographic/Co-ordinates instead of pixel Coords?

Question

How do I get the 1D profile from an Array given geographic coordinates instead of pixel coordinates. I have radiance data in Geo-coordinates (NASA GOLD Mission). The data for each hemisphere is stored in a separate file. I want to take a 1D profile of the spliced hemispheres along a curve (10 degNorth dip latitude) which is in Geo coords,

I tried to adapt the answer here using map_coordinates but I realised it is in pixel coordinates that is why I believe the profile does not correspond to the map; To get the 1D profile I plotted the output from mapcoords against the input curve (X) values so as to get the variation against longitude. The other profile in yellow also does not correspond to the data in the map; When I plot the data using imshow it does not look the same as on the map;

The data and Python Code is here Data+Code;

from netCDF4 import Dataset  
import glob, matplotlib,os;matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy.ma as ma

Mag = pd.read_csv(r"\mag_cords.csv")
X_dip10s,Y_dip10s = Mag.LON10S.to_numpy(),Mag.LAT10S.to_numpy()
g = Dataset('GOLD_L1C_CH*23*.nc','r'); 
wavelength = g.variables['WAVELENGTH'][:]
radiance = g.variables['RADIANCE'][:]
lat = g.variables['REFERENCE_POINT_LAT'][:] 
lon = g.variables['REFERENCE_POINT_LON'][:] 
g.close()
O5s_ids = np.argwhere((134.3 <= wavelength[50,25,:]) & (wavelength[50,25,:] <= 137.7))
O5s = np.array(np.nansum(radiance[:,:,O5s_ids],axis = 2))*.04 # integrate under the peak!
O5s = np.transpose(O5s[:,:,0])
x=lon[9:-8,9:-8];y=lat[9:-8,9:-8];z=O5s[9:-8,9:-8].T;z=abs(z)#To remove Nans
z2=np.ma.masked_array(z, (z > 30)) 

Zs=scipy.ndimage.map_coordinates(np.transpose(z),
np.vstack((X_dip10s,Y_dip10s)),
mode="nearest")
Xs=scipy.ndimage.map_coordinates(np.transpose(x),
np.vstack((X_dip10s,Y_dip10s)),mode="nearest") 
                                                                                                                         
fig=plt.figure(figsize=(10,10))
axarr = fig.add_subplot(211,projection=ccrs.PlateCarree())
ax2 = fig.add_subplot(212)
grid_lines =axarr.gridlines(draw_labels=True'); axarr.coastlines()
cs1 = axarr.contourf(x,y, z2, transform=ccrs.PlateCarree(),cmap='jet',levels=50)
ax2.plot(X_dip10s,Zs, 'k-', Xs ,Zs, 'k-')


  [1]: https://i.stack.imgur.com/66C4P.png

Answer 1

interesting question!

I'm the dev of EOmaps, a library intended to simplify creating exactly this kind of custom geo-data analysis widgets, so I got interested :-)

I had a quick look at your data... its a quite strange netcdf and I have no clue what defines the line you use to identify the points so I just re-used your code to extract the data.

Here's a quick solution that uses a KDTree to identify the closest points with respect to the line-points... Its not 100% exact since it identifies those points with respect to a max. distance defined in degrees and not actual distance on the globe, but it works well enough and I think you should be able to take it from here if you need something more exact...

Its also important to notice that this method will show you the actual data and not some interpolated points created by the algorithm used to display the data! (so the line will not be nice and smooth but it will only show the actual data values closest to the line)

from pathlib import Path
from itertools import chain

from scipy.spatial import KDTree
import xarray as xar
import pandas as pd
import numpy as np

from eomaps import Maps


# --------------------------------------------- read the data
p = Path(r"GOLD_data-20230319T161539Z-001\GOLD_data")

with xar.open_dataset(p / "GOLD_L1C_CHA_NI1_2018_286_23_10_v04_r01_c01.nc") as ncfile:
    lon, lat = ncfile.REFERENCE_POINT_LON.values, ncfile.REFERENCE_POINT_LAT.values
    wavelength = ncfile.WAVELENGTH.values
    radiance = ncfile.RADIANCE.values

O5s_ids = np.argwhere((134.3 <= wavelength[50,25,:]) & (wavelength[50,25,:] <= 137.7))
O5s = np.array(np.nansum(radiance[:,:,O5s_ids],axis = 2))*.04 # integrate under the peak!
O5s = np.transpose(O5s[:,:,0])
x=lon[9:-8,9:-8]
y=lat[9:-8,9:-8]
z=O5s[9:-8,9:-8].T


Mag = pd.read_csv(p / "mag_cords.csv")
X_dip10s,Y_dip10s = Mag.LON10S.to_numpy(),Mag.LAT10S.to_numpy()
line = np.column_stack((X_dip10s, Y_dip10s))


# --------------------------------------------- create a new map
m = Maps(ax=211, crs=4326)
m.set_extent((-61.8,-20.5,-1.5,24.5))
m.add_feature.preset.coastline()
m.add_feature.preset.ocean()

# plot the data
m.set_data(data=z, x=x, y=y, crs=4326)
m.set_shape.ellipses(radius=0.35)
m.set_classify.UserDefined(bins=np.linspace(0, 20, 11))
m.plot_map(cmap="viridis", vmin=0, vmax=28)
m.add_colorbar(orientation="vertical")

# add a normal matplotlib axis (and share x-limits)
ax = m.f.add_subplot(212)
ax.sharex(m.ax)
# turn off navigation on the normal matplotlib plot
ax.set_navigate(False)

# position axes 
m.apply_layout(
    {'figsize': [6.95, 4.8],
     '0_map': [0.1, 0.44678, 0.5875, 0.53552],
     '1_cb': [0.75, 0.04546, 0.2125, 0.95484],
     '1_cb_histogram_size': 0.8,
     '2_': [0.1, 0.10859, 0.5875, 0.32578]}
    )

# setup a new map-layer (will be used to draw identified points)
m2 = m.new_layer()

# --------------------------------------------- identify points along line
tree = KDTree(np.column_stack((x.ravel(), y.ravel())))

def get_lineoffset(pos, **kwargs):
    # get the offset required to shift the line to the click position
    x, y = line.T
    minxid = np.argmin(abs(pos[0] - x))
    offset = pos[1] - y[minxid]
    return offset

def identify_points(pts, max_dist):
    # identify closest points with respect to the line
    tree2 = KDTree(pts)
    q = tree2.query_ball_tree(tree, max_dist)
    q = list(chain(*[i for i in q if len(i) > 0]))
    coords = tree.data[q].T
    data = z.ravel()[q]
    return coords, data

def indicate_found_points(pts, offset=0, max_dist=1):
    # offset line
    pts = pts + [0, offset]
    # identify closest points
    coords, data = identify_points(pts, max_dist=max_dist)
    
    # plot identified points on the map (replace on each click)
    m2.set_data(data, coords[0], coords[1])
    m2.set_shape.ellipses(radius=0.2)
    m2.plot_map(set_extent=False, zorder=2, fc="r")
    
    # sort found points so we can plot a line
    sortp = np.argsort(coords[0])
    lx, ly = coords[0][sortp], data[sortp]    
    
    # plot a line on the map
    l, = m.ax.plot(*pts.T, c="r", lw=0.5)
    # remove it on next pick
    m.cb.pick.add_temporary_artist(l)
    
    # plot data on normal matpltoib plot
    l2, = ax.plot(lx, ly, lw=0.25, c="k", zorder=0)
    # color points same as in the map
    sc = ax.scatter(lx, ly, c=ly, s=10, norm=m._norm, cmap=m._cbcmap)
    # remove it on next pick
    m.cb.pick.add_temporary_artist(l2)
    m.cb.pick.add_temporary_artist(sc)

    ax.relim()
    ax.autoscale_view()
    m.redraw()

def callback(pos, max_dist=1, **kwargs):
    offset = get_lineoffset(pos)
    indicate_found_points(line, offset=offset, max_dist=max_dist)

# attach the callback to the map
m.cb.pick.attach(callback, max_dist=0.75)