计算BaseMap创建的轮廓线内的非项目

Question

我目前正在尝试使用BaseMap在Mollweide地图投影上确定Specfic轮廓线中的区域。具体而言，我正在寻找的是天文学上的天文学事件的各种可靠间隔的区域。该图如下所示：

幸运的是，a 类似的问题< /a>在此之前已经得到了回答。答案中概述的方法也能够说明轮廓内的孔，这对于我的用例来说是必要的。下面提供了我针对此特定方法的改编代码：

# 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)

据我所知，此方法可很好地工作...除了一个小皱纹：这可以使用投影坐标单元计算该区域。我不完全确定在这种情况下这些单位是什么，但是它们绝对不是学位 ² （计算的区域处于10 ¹³ 单位 ² ...也许单位是像素？）。如前所述，我正在寻找的是如何计算全局坐标单元中的等效区域，即以度为单位 ² 。

是否有一种方法可以将预计域中计算的区域转换为平方度的全球域？也许有一种方法可以修改此方法，以便从一开始就确定该区域 ² ？

将不胜感激！

Answer 1

对于遇到这个问题的任何人，虽然我没有找到一种直接将投影区域转换回全球域的方法，但我确实通过转换轮廓路径顶点开发了一个新的解决方案（但是这次是在LAT/中定义的LON坐标系）通过保留正弦投影的区域：

其中φ是纬度，λ是经度，而λ ₀ 是中央子午线的经度。

此扁平投影意味着您只能将包装齐头发地确定由投影顶点定义的多边形区域（以1个单位的半径为1个或更简单的steradians）。将该数字乘以（180/π） ² 将为您的轮廓以平方度为您的面积。

幸运的是，仅需要对OP中提到的代码进行少量调整才能实现这一目标。最终代码如下：

# 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)