1

我生成了一个 Voronoi 图并计算了面积，然后计算了形成的多边形的 1/面积（通量）并找到了平均通量。我现在想按通量对我的多边形进行颜色编码（我只想用通量>平均通量对多边形着色）

#Voronoi algorithm

list_ra=tbdata[cut3]['RA']

list_dec=tbdata[cut3]['DEC']

points=np.column_stack((list_ra, list_dec))

vor=Voronoi(points)

#print vor.vertices

#print vor.regions

#print vor.ridge_vertices
#print vor.ridge_points

print len(vor.point_region)

plt.plot(points[:, 0], points[:, 1], '.')

plt.plot(vor.vertices[:, 0], vor.vertices[:, 1], '.', markersize=0.25)

plt.xlim(150.02, 150.21); plt.ylim(2.26, 2.45)

for simplex in vor.ridge_vertices:

     simplex = np.asarray(simplex)

     if np.all(simplex >= 0):

      plt.plot(vor.vertices[simplex, 0], vor.vertices[simplex, 1], 'k-')

center = points.mean(axis=0)

for pointidx, simplex in zip(vor.ridge_points, vor.ridge_vertices):

     simplex = np.asarray(simplex)

     if np.any(simplex < 0):

         i = simplex[simplex >= 0][0] # finite end Voronoi vertex

         t = points[pointidx[1]] - points[pointidx[0]]  # tangent

         t = t / np.linalg.norm(t)

         n = np.array([-t[1], t[0]]) # normal

         midpoint = points[pointidx].mean(axis=0)

         far_point = vor.vertices[i] + np.sign(np.dot(midpoint - center, n)) * n * 100

         plt.plot([vor.vertices[i,0], far_point[0]],

                  [vor.vertices[i,1], far_point[1]], 'k--')

plt.xlabel(r' RA/deg ',size=12)

plt.ylabel(r' DEC/deg ',size=12)

plt.savefig('Voronoi_10x10', dpi=300)

plt.show()

###############################################################################################
#Polygon areas

total=0

totals=[]

for i in range (len(vor.point_region)):

    polygon_coordinates=[vor.vertices[y] for y in vor.regions[vor.point_region[i]]]

    list_polygon_coordinates=np.asarray(polygon_coordinates)

    r=list_polygon_coordinates[:, 0]

    d=list_polygon_coordinates[:, 1]

    def PolyArea(r,d): 
        return 0.5*np.abs(np.dot(r,np.roll(d,1))-np.dot(d,np.roll(r,1)))

    total=PolyArea(r,d)

    totals.append(total)

totaltotal=0

for i in range(len(vor.point_region)):

    totaltotal+=totals[i]

print ('totaltotal', totaltotal)

average_area=totaltotal/1533

print ('average_area', average_area)

################################################################################################
#Flux distribution

area=np.array(totals)

print ('area', area)

flux=1/area

print ('flux', flux)

totaltotal_flux=0

for i in range(1533):

    totaltotal_flux+=flux[i]

print ('totaltotal_flux', totaltotal_flux)

average_flux=totaltotal_flux/1533

print ('average_flux', average_flux)
4

0 回答 0