pythonpython-3.xpandasmappingpygmt

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


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:

Should be an elevation map of the Lunar South Pole (All Longitudes, -90 - -80 Latitude)

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?


Solution

  • 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