Calculate the non-projected area inside a contour line created by Basemap

I am currently trying to determine the area inside specfic contour lines on a Mollweide map projection using Basemap. Specifically, what I'm looking for is the area of various credible intervals in square degrees (or degrees²) of an astronomical event on the celestial sphere. The plot is shown below:

Fortunately, a similar question has already been answered before that helps considerably. The method outlined in the answer is able to account for holes within the contour as well which is a necessity for my use case. My adapted code for this particular method is provided below:

# generate a regular lat/lon grid.
nlats = 300; nlons = 300; delta_lon = 2.*np.pi/(nlons-1); delta_lat = np.pi/(nlats-1)
lats = (0.5*np.pi-delta_lat*np.indices((nlats,nlons))[0,:,:])
lons = (delta_lon*np.indices((nlats,nlons))[1,:,:] - np.pi)

map = Basemap(projection='moll',lon_0=0, celestial=True)

# compute native map projection coordinates of lat/lon grid
x, y = map(lons*180./np.pi, lats*180./np.pi)    

areas = []
cred_ints = [0.5,0.9]

for k in range(len(cred_ints)):

    cs = map.contourf(x,y,p1,levels=[0.0,cred_ints[k]]) ## p1 is the cumulative distribution across all points in the sky (usually determined via KDE on the data)
    
    ##organizing paths and computing individual areas
    paths = cs.collections[0].get_paths()
    #help(paths[0])
    area_of_individual_polygons = []
    for p in paths:
        sign = 1  ##<-- assures that area of first(outer) polygon will be summed
        verts = p.vertices
        codes = p.codes
        idx = np.where(codes==Path.MOVETO)[0]
        vert_segs = np.split(verts,idx)[1:]
        code_segs = np.split(codes,idx)[1:]
        for code, vert in zip(code_segs,vert_segs):

            ##computing the area of the polygon
            area_of_individual_polygons.append(sign*Polygon(vert[:-1]).area)
            sign = -1 ##<-- assures that the other (inner) polygons will be subtracted

    ##computing total area        
    total_area = np.sum(area_of_individual_polygons)
    
    print(total_area)
    
    areas.append(total_area)

print(areas)

As far as I can tell this method works beautifully... except for one small wrinkle: this calculates the area using the projected coordinate units. I'm not entirely sure what the units are in this case but they are definitely not degrees² (the calculated areas are on the order of 10¹³ units²... maybe the units are pixels?). As alluded to earlier, what I'm looking for is how to calculate the equivalent area in the global coordinate units, i.e. in degrees².

Is there a way to convert the area calculated in the projected domain back into the global domain in square degrees? Or perhaps is there a way to modify this method so that it determines the area in degrees² from the get go?

Any help will be greatly appreciated!

For anyone that comes across this question, while I didn't figure out a way to directly convert the projected area back into the global domain, I did develop a new solution by transforming the contour path vertices (but this time defined in the lat/lon coordinate system) via an area preserving sinusoidal projection:

where φ is the latitude, λ is the longitude, and λ₀ is the longitude of the central meridian.

This flat projection means you can just use the package Shapely to determine the area of the polygon defined by the projected vertices (in square units for a radius of 1 unit, or more simply steradians). Multiplying this number by (180/π)² will give you the area in square degrees for the contour in question.

Fortunately, only minor adjustments to the code mentioned in the OP was needed to achieve this. The final code is provided below:

# generate a regular lat/lon grid.
nlats = 300; nlons = 300; 
delta_lat = np.pi/(nlats-1); delta_lon = 2.*np.pi/(nlons-1);    
lats = (0.5*np.pi-delta_lat*np.indices((nlats,nlons))[0,:,:])
lons = (delta_lon*np.indices((nlats,nlons))[1,:,:])    

### FOLLOWING CODE DETERMINES CREDIBLE INTERVAL SKY AREA IN DEG^2  ###

# collect and organize contour data for each credible interval

cred_ints = [0.5,0.9]    
ci_areas = []

for k in range(len(cred_ints)):        

    cs = plt.contourf(lons,lats,p1,levels=[0,cred_ints[k]])   ## p1 is the cumulative distribution across all points in the sky (usually determined via KDE of the dataset in question)
    paths = cs.collections[0].get_paths()
    
    ##organizing paths and computing individual areas

    area_of_individual_polygons = []

    for p in paths:
        
        sign = 1  ##<-- assures that area of first(outer) polygon will be summed
        vertices = p.vertices
        codes = p.codes
        idx = np.where(codes==Path.MOVETO)[0]    
        verts_segs = np.split(vertices,idx)[1:]            
        
        for verts in verts_segs:                
            
            # transforming the coordinates via an area preserving sinusoidal projection 

            x = (verts[:,0] - (0)*np.ones_like(verts[:,0]))*np.cos(verts[:,1])
            y = verts[:,1]

            verts_proj = np.stack((x,y), axis=1)        

            ##computing the area of the polygon                  
            area_of_individual_polygons.append(sign*Polygon(verts_proj[:-1]).area)                
            sign = -1 ##<-- assures that the other(inner) polygons/holes will be subtracted

    ##computing total area        
    total_area = ((180/np.pi)**2)*np.sum(area_of_individual_polygons)
    
    ci_areas.append(total_area)