What am I doing wrong when trying to plot an elevation map in pyGMT using Pandas?

Question

I have a data file for elevation data from the Japanese Lunar Orbiter SELENE (nicknamed Kaguya). The data is found here: https://data.darts.isas.jaxa.jp/pub/pds3/sln-l-lalt-5-topo-gt-sp-num-v2.0/data/LALT_GT_SP_NUM.TAB

It is a 450 MB data of Long Lat Alt values running from -79ish latitude to -90ish latitude. and 0 360 longitude.

I managed to clean up the data (remove extra spaces) and put it into a pandas data frame as well as plot it using pyGMT, which is a python module for interfacing with the Generic Mapping Tools program.

My plot however is messed up, it looks like this:

This came out incorrectly, and near the center there are two empty gray regions that should be populated with elevation data. I will include my code below for plotting the data. The data can be accessed fairly easily but its a large file so I will put in a few example rows of it as well for anyone that does not want to download the file themselves.

"df" is the pandas dataframe I pass to the plot function.

def plotGrid(df):
    startTimeOA = time.time()
    grid = pygmt.xyz2grd(
        data=df,
        projection="E0/-90/20c",
        region='0/360/-90/-80',
        spacing=[0.03125, 0.0078125]
    )
    print(f'Took {time.time() - startTimeOA:06.3f} seconds to create grid from relevant Pandas dataFrame.')
    startTime = time.time()
    fig = pygmt.Figure()
    fig.grdimage(grid=grid)
    fig.show()
    print(f'Took {time.time() - startTimeOA:06.3f} seconds overall to plot Lunar South Pole Map '
          f'from relevant LALT data.')
    print('Plotted Long Lat Alt elevation data in pyGMT for relevant LALT data.\n')
    return fig

Here is an example of what the data itself looks like. The columns are long, lat, and altitude respectively in degrees and kilometers:

0.015625 -79.00390625 2.114
0.046875 -79.00390625 2.123
0.078125 -79.00390625 2.133
0.109375 -79.00390625 2.144
0.140625 -79.00390625 2.160
0.171875 -79.00390625 2.177
0.203125 -79.00390625 2.195
0.234375 -79.00390625 2.213
0.265625 -79.00390625 2.230

For "projection", I have my map projection as "E" which is Azimuthal Equidistant. It's centered at 0 longitude and -90 latitude.

The 20c after that specifies the map size as 20 centimeters I believe.

My questions are, why does the map look like this if the data is equally spaced in the file?

Longitude varies in increments of 0.03125 degrees and latitude in increments of 0.0078125 degrees. I have set the spacing to that as well.

So my questions are: What am I doing wrong here in trying to create an elevation map from this data? What do I need to do to get rid of these weird radial gaps along the top and bottom edges as well as the gaps at the top and bottom of the center of the map?

Answer 1

I think you have to set the registration parameter of xyz2grd to "p" for pixel registration.

import pandas as pd
import pygmt as gmt

# %%
# >>> huge file - takes some time <<<
url_data = "https://data.darts.isas.jaxa.jp/pub/pds3/sln-l-lalt-5-topo-gt-sp-num-v2.0/data/LALT_GT_SP_NUM.TAB"
data = pd.read_csv(url_data, sep="\s+", header=None, names=["lon", "lat", "alt"])

# %%
region = [0, 360, -90, -80]
registration = "p"

grid = gmt.xyz2grd(
    data=data,
    region=region,
    spacing=[0.03125, 0.0078125],
    registration=registration,
)

fig = gmt.Figure()
fig.grdimage(
    grid=grid,
    projection="E0/-90/20c",
    region=region,
    frame=True,
    cmap="batlow",
)
fig.show()

# fig.savefig(fname=f"map_moon_{registration}.png")

Output figures:

pixel registration	gridline registration

Answer 2

I'm unfortunately not familiar with pyGMT, but your data looked really interesting so I thought I'll give it a quick try and see if EOmaps can do it (I'm the dev :-)

I downloaded the data and read it like this:

import pandas as pd
data = pd.read_csv("LALT_GT_SP_NUM.TAB", sep="\s+", header=None, names=["lon", "lat", "alt"])
# Convert longitude form 0/360 to -180/180
data.lon = (data.lon + 180)%360 -180

And here's an example how you can quickly visualize your data:

from eomaps import Maps
m = Maps(crs=Maps.CRS.AzimuthalEquidistant(central_latitude=-90, central_longitude=0))
m.set_data(data)
m.set_shape.shade_points()
m.plot_map()
m.add_colorbar(label="Japanese Lunar Orbiter SELENE Altitude")

(I used "shade_points" as a plot-shape to avoid problems with the data-size of ~50M datapoints and to avoid having to bother converting it to 2D prior to plotting... even though that might be a good idea on the long run...)

... and the nice thing is that I've implemented access to the OpenPlanetary WebMap layers sometime ago so you can directly compare your data to their Moon basemap!

m.add_wms.OpenPlanetary.Moon.add_layer.basemap_layer(layer="OP_basemap")
m.cb.click.attach.peek_layer("OP_basemap")