How to plot data on a basemap using matplotlib basemap

Question

Two sections of my code are giving me trouble, I am trying to get the basemap created in this first section here:

#Basemap
epsg = 6060; width = 2000.e3; height = 2000.e3 #epsg 3413. 6062
m=Basemap(epsg=epsg,resolution='l',width=width,height=height) #lat_ts=(90.+35.)/2.
m.drawcoastlines(color='white')
m.drawmapboundary(fill_color='#99ffff')
m.fillcontinents(color='#cc9966',lake_color='#99ffff')
m.drawparallels(np.arange(10,70,20),labels=[1,1,0,0])
m.drawmeridians(np.arange(-100,0,20),labels=[0,0,0,1])
plt.title('ICESAT2 Tracks in Greenland')
plt.figure(figsize=(20,10))

Then my next section is meant to plot the data its getting from a file, and plot these tracks on top of the Basemap. Instead, it creates a new plot entirely. I have tried rewording the secondary plt.scatter to match Basemap, such as m.scatter, m.plt, etc. But it only returns with “RuntimeError: Can not put single artist in more than one figure” when I do so. Any ideas on how to get this next section of code onto the basemap? Here is the next section, focus on the end to see where it is plotting.

icesat2_data[track] = dict() # creates a sub-dictionary, track
    icesat2_data[track][year+month+day] = dict() # and one layer more for the date under the whole icesat2_data dictionary
    icesat2_data[track][year+month+day] = dict.fromkeys(lasers)
    for laser in lasers: # for loop, access all the gt1l, 2l, 3l
        if laser in f:    
            lat = f[laser]["land_ice_segments"]["latitude"][:] # data for a particular laser's latitude.
            lon = f[laser]["land_ice_segments"]["longitude"][:] #data for a lasers longitude
            height = f[laser]["land_ice_segments"]["h_li"][:] # data for a lasers height
            quality = f[laser]["land_ice_segments"]["atl06_quality_summary"][:].astype('int')
        
                # Quality filter
            idx1 = quality == 0 # data dictionary to see what quality summary is
            #print('idx1', idx1)
        
                # Spatial filter
            idx2 = np.logical_and( np.logical_and(lat>=lat_min, lat<=lat_max), np.logical_and(lon>=lon_min, lon<=lon_max) )
        
            idx = np.where( np.logical_and(idx1, idx2) ) # combines index 1 and 2 from data quality filter. make sure not empty. if empty all data failed test (low quality or outside box)

            icesat2_data[track][year+month+day][laser] = dict.fromkeys(['lat','lon','height']) #store data, creates empty dictionary of lists lat, lon, hi, those strings are the keys to the dict.
            icesat2_data[track][year+month+day][laser]['lat'] = lat[idx] # grabbing only latitudes using that index of points with good data quality and within bounding box
            icesat2_data[track][year+month+day][laser]['lon'] = lon[idx] 
            icesat2_data[track][year+month+day][laser]['height'] = height[idx]
         
            if lat[idx].any() == True and lon[idx].any() == True:
                x, y = transformer.transform(icesat2_data[track][year+month+day][laser]['lon'], \
                                        icesat2_data[track][year+month+day][laser]['lat']) 

                plt.scatter(x, y, marker='o', color='#000000')

Currently, they output separately, like this:

Answer 1

Not sure if you're still working on this, but here's a quick example I put together that you might be able to work with (obviously I don't have the data you're working with). A couple things that might not be self-explanatory - I used m() to transform the coordinates to map coordinates. This is Basemap's built-in transformation method so you don't have to use PyProj. Also, setting a zorder in the scatter function ensures that your points are plotted above the countries layer and don't get hidden underneath.

#Basemap

epsg = 6060; width = 2000.e3; height = 2000.e3 #epsg 3413. 6062
plt.figure(figsize=(20,10))
m=Basemap(epsg=epsg,resolution='l',width=width,height=height) #lat_ts=(90.+35.)/2.
m.drawcoastlines(color='white')
m.drawmapboundary(fill_color='#99ffff')
m.fillcontinents(color='#cc9966',lake_color='#99ffff')
m.drawparallels(np.arange(10,70,20),labels=[1,1,0,0])
m.drawmeridians(np.arange(-100,0,20),labels=[0,0,0,1])
plt.title('ICESAT2 Tracks in Greenland')


for coord in [[68,-39],[70,-39]]:
    lat = coord[0]
    lon = coord[1]
    x, y = m(lon,lat)
    m.scatter(x,y,color='red',s=100,zorder=10)
              
plt.show()

Answer 2

I think you might need:

plt.figure(figsize(20,10))

before creating the basemap, not after. As it stands it's creating a map and then creating a new figure after that which is why you're getting two figures.

Then your plotting line should be m.scatter() as you mentioned you tried before.