我在写一个代码,比较天文图上像素的通量和另一张地图上相应区域的通量。两张地图都是大量的数据。你知道吗
为了做到这一点,我需要将第一张地图上的像素索引(Av)转换为它们在天空坐标上的等效索引,然后将这些天空坐标转换为它们在第二张地图上的像素索引(CO)。然后,我缩放第二个贴图的通量以匹配第一个贴图的值。在那之后,我必须继续处理数据。你知道吗
问题是,由于第一张地图上有数千个像素,代码要花很长时间才能完成它应该做的事情,这对于故障排除来说是一个麻烦。我发现这部分代码中最慢的是for循环。你知道吗
有没有什么方法可以遍历numpy数组,能够处理索引并计算每个像素的数据,比for循环更快?有没有更好的方法?你知道吗
在伪代码中,我的代码是这样的:
for pixel i,j in 1st map:
sky_x1,sky_y1 = pixel_2_skycoord(i,j)
i2,j2 = skycoord_2_pixel(sky_x1,sky_y1)
Avmap.append(Avflux[i,j])
COmap.append(COflux[i2,j2]*scale)
实际代码为:
for i in xrange(0,sAv_y-1):
for j in xrange(0,sAv_x-1):
if not np.isnan(Avdata[i,j]):
y,x=wcs.utils.skycoord_to_pixel(wcs.utils.pixel_to_skycoord(i,j,wAv,0),wcs=wCO)
x=x.astype(int)+0 #the zero is because i don't understand the problem with numpy but it fixes it anyway
y=y.astype(int)+0 #i couldn't get the number from an array with 1 value but adding zero resolves it somehow
COflux=COdata[x,y]
ylist.append(Avdata[i,j])
xlist.append(COflux*(AvArea/COArea))
罪魁祸首是两个for循环。Numpy有许多函数可以防止使用for循环来允许快速编译代码。诀窍是将代码矢量化。你知道吗
您可以查看numpy的
meshgrid
函数,将此数据转换为向量化形式,然后可以使用this SO question之类的函数对该向量应用任意函数。你知道吗大致如下:
您必须避免循环,并在底层的C代码中、在Numpy或Astropy中进行繁重的计算,以便进行天空/像素转换。有几个选项可以用
astropy.wcs
来实现这一点。你知道吗第一个是用
SkyCoord
。我们首先为像素索引创建一个值网格:现在我们可以从像素索引创建
SkyCoord
对象(这是Numpy数组子类),并使用wcs:注意,这是使用
wcs.utils.skycoord_to_pixel
。这个对象还有一个方法,可以用wcs投影到像素。出于实际目的,我在这里也将这样做:我们得到一个新的x和y索引的浮点值元组。因此,您需要对这些值进行舍入,并将其转换为int以用作数组索引。你知道吗
第二个选项是使用较低级别的函数,例如
wcs.pixel_to_world_values
和wcs.world_to_pixel_values
,它接受Nx2数组并返回:相关问题 更多 >
编程相关推荐