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

Question

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!

Answer 1

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 modifying the original approach and utilizing Green's theorem with the contour path information. First, some background information (the code is provided at the end of the answer for anyone not interested):

The solid angle encompassed by a region R on a sphere (credible intervals on the celestial sphere in this context) is given by the following:

Using Green's theorem, i.e.

we can convert our surface integral into a contour integral over the region's boundary C. We just need to solve for the functions Q and P which is relatively straightforward,

Without loss of generality, the simplest solution to this equation is

Plugging these functions into Green's theorem and equating it with the total solid angle gives us our final equation:

Using the contour path information from matplotlib.contourf, and the sign specifying if a path is a hole or not, we can now easily determine the total solid angle any particular contour C encompasses. The only complication is accounting for negative areas (in this case contours that traverse the lower hemisphere). However, this is easily resolved by splitting up the path information for a contour into separate paths for each hemisphere and negating the areas obtained for the paths in the lower hemisphere.

The code that actually carries out the calculation 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,:,:])


### 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]    
all_lines = []

for k in range(len(cred_ints)):       

    cs = plt.contourf(lons,lats,p1,levels=[0.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()

    lines = []

    for p in paths:            

        sign = 1  ##<-- assures that area of first(outer) paths 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:]            

        temp = []

        for code, vert in zip(code_segs,vert_segs):                

            temp.append((sign,vert))                
            sign = -1 ##<-- assures that the other (inner) paths/holes will be subtracted

        lines.append(temp)

    all_lines.append(lines) 


# use contour data to calculate relevant areas

ci_areas = []

for k in range(len(cred_ints)):  ### number of credible intervals 

    ci_contours = all_lines[k]

    contour_areas = []

    for j in range(len(ci_contours)):  ### number of different contours for current credible interval       

        contour_paths = ci_contours[j]

        path_areas = [] 

        for i in range(len(contour_paths)):  ### number of unique paths that make up current contour (includes holes)

            areas = []

            sign = contour_paths[i][0]
            path = contour_paths[i][1]

            pos_idx = np.where(path[:,1]>=0)[0] ## upper hemisphere paths
            pos_path = np.vstack(np.split(path[pos_idx],np.where(np.diff(pos_idx)!=1)[0]+1))

            neg_idx = np.where(path[:,1]<0)[0]  ## lower hemisphere paths
            temp_neg_path = np.vstack(np.split(path[neg_idx],np.where(np.diff(neg_idx)!=1)[0]+1))                        
            neg_path = np.roll(temp_neg_path, (np.shape(temp_neg_path)[0]//4)*2)  ##cycle array so the cutoff occurs near the middle and not at the beginning/end (this leads to spurious values for some reason)


            ## the np.abs step in the following lines should be redundant but I've kept it in for peace of mind :)   

            posA = sign*np.abs(np.sum(-np.cos(pos_path[:,1])[:-1]*np.diff(pos_path[:,0])))

            areas.append(posA)

            negA = sign*np.abs(np.sum(np.cos(neg_path[:,1])[:-1]*np.diff(neg_path[:,0])))

            areas.append(negA)

            path_areas.append(np.sum(areas))

        single_contour_area = np.sum(path_areas)

        contour_areas.append(single_contour_area)

    total_area = np.sum(contour_areas)

    ci_areas.append(total_area*(180/np.pi)**2)

print(ci_areas)