计算底图创建的等高线内的非投影区域
Calculate the non-projected area inside a contour line created by Basemap
我目前正在尝试使用底图确定 Mollweide 地图投影上特定等高线内的区域。具体来说，我正在寻找的是天球上天文事件的平方度（或度数2）的各种可信区间的面积。剧情如下图：
幸运的是，similar question 之前已经有人回答了，这很有帮助。答案中概述的方法也能够解决轮廓内的孔洞，这对我的用例来说是必需的。下面提供了我为此特定方法改编的代码：
# 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)
据我所知，这种方法效果很好……除了一个小问题：它使用投影坐标单位计算面积。我不完全确定这种情况下的单位是什么，但它们绝对不是度数2（计算出的面积约为 1013 单位2...也许单位是像素？）。如前所述，我正在寻找的是如何计算全局坐标单位中的等效面积，即度数2.
有没有办法将投影域中计算的面积转换回全局域的平方度？或者也许有一种方法可以修改此方法，以便它从一开始就以度数 2 为单位确定面积？
任何帮助将不胜感激！
对于遇到这个问题的任何人，虽然我没有找到将投影区域直接转换回全局域的方法，但我确实通过修改原始方法并利用格林定理开发了一个新的解决方案轮廓路径信息。首先，一些背景信息（代码在答案末尾提供给任何不感兴趣的人）：
球体上的区域 R 所包含的立体角（在此上下文中天球上的可信间隔）由以下给出：
利用格林定理，即
我们可以将我们的表面积分转换为区域边界C上的轮廓积分。我们只需要求解函数Q和P，这相对简单，
不失一般性，这个方程最简单的解是
将这些函数代入格林定理，并将其与总立体角等同，得出最终方程：
使用来自 matplotlib.contourf 的轮廓路径信息，以及指定路径是否为孔的符号，我们现在可以轻松确定任何特定轮廓 C 包含的总立体角。唯一复杂的是计算负区域（在本例中为横跨下半球的等高线）。然而，这很容易解决，方法是将等高线的路径信息拆分为每个半球的单独路径，并否定为下半球路径获得的区域。
实际执行计算的代码如下：
# 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)
我目前正在尝试使用底图确定 Mollweide 地图投影上特定等高线内的区域。具体来说，我正在寻找的是天球上天文事件的平方度（或度数2）的各种可信区间的面积。剧情如下图：
幸运的是，similar question 之前已经有人回答了，这很有帮助。答案中概述的方法也能够解决轮廓内的孔洞，这对我的用例来说是必需的。下面提供了我为此特定方法改编的代码：
# 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)
据我所知，这种方法效果很好……除了一个小问题：它使用投影坐标单位计算面积。我不完全确定这种情况下的单位是什么，但它们绝对不是度数2（计算出的面积约为 1013 单位2...也许单位是像素？）。如前所述，我正在寻找的是如何计算全局坐标单位中的等效面积，即度数2.
有没有办法将投影域中计算的面积转换回全局域的平方度？或者也许有一种方法可以修改此方法，以便它从一开始就以度数 2 为单位确定面积？
任何帮助将不胜感激！
对于遇到这个问题的任何人，虽然我没有找到将投影区域直接转换回全局域的方法，但我确实通过修改原始方法并利用格林定理开发了一个新的解决方案轮廓路径信息。首先，一些背景信息（代码在答案末尾提供给任何不感兴趣的人）：
球体上的区域 R 所包含的立体角（在此上下文中天球上的可信间隔）由以下给出：
利用格林定理，即
我们可以将我们的表面积分转换为区域边界C上的轮廓积分。我们只需要求解函数Q和P，这相对简单，
不失一般性，这个方程最简单的解是
将这些函数代入格林定理，并将其与总立体角等同，得出最终方程：
使用来自 matplotlib.contourf 的轮廓路径信息，以及指定路径是否为孔的符号，我们现在可以轻松确定任何特定轮廓 C 包含的总立体角。唯一复杂的是计算负区域（在本例中为横跨下半球的等高线）。然而，这很容易解决，方法是将等高线的路径信息拆分为每个半球的单独路径，并否定为下半球路径获得的区域。
实际执行计算的代码如下：
# 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)