-
Notifications
You must be signed in to change notification settings - Fork 8
/
Copy pathstructurefunction.py
232 lines (183 loc) · 8.36 KB
/
structurefunction.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
# This file is part of MAGNETAR, the set of magnetic field analysis tools
#
# Copyright (C) 2013-2019 Juan Diego Soler
import sys
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.convolution import convolve, convolve_fft
from astropy.convolution import Gaussian2DKernel
from astropy.stats.circstats import *
from nose.tools import assert_equal, assert_true
from scipy.ndimage import gaussian_filter
from bvisual import *
from poltools import *
from astropy.wcs import WCS
def s2(Qmap, Umap, lags=[4.0], s_lag=1.0, mask=None, header=None):
assert Qmap.shape == Umap.shape, "Dimensions of Qmap and Umap must match"
sz=np.shape(Qmap)
QmapRAND=(np.random.rand(sz[0],sz[1])-0.5)
UmapRAND=(np.random.rand(sz[0],sz[1])-0.5)
if mask is None:
mask=np.ones_like(Qmap)
assert Qmap.shape == mask.shape, "Dimensions of mask and Qmap must match"
Qmap[(mask==0.).nonzero()]=np.nan
Umap[(mask==0.).nonzero()]=np.nan
# ------------------------------------------------------------------------------------------------
if header is None:
posx, posy=np.meshgrid(np.arange(0,sz[0]), np.arange(0,sz[1]))
else:
ra=header['CRVAL1']+header['CDELT1']*(np.arange(header['NAXIS1'])-header['CRPIX1'])
dec=header['CRVAL2']+header['CDELT2']*(np.arange(header['NAXIS2'])-header['CRPIX2'])
posx, posy=np.meshgrid(ra, dec)
print('Calculation of x positions')
x1, x2 =np.meshgrid(posx.ravel(), posx.ravel())
deltax=(x1-x2)
x1=None; x2=None
print('Calculation of y positions')
y1, y2 =np.meshgrid(posy.ravel(), posy.ravel())
deltay=(y1-y2)
y1=None; y2=None
print('Calculation of distances')
dist=np.sqrt(deltax**2+deltay**2)
deltax=None; deltay=None;
# ----------------------------------------------------------------------------------------------
stwo=np.nan
stwoRAND=np.nan
stwoarr=np.nan*np.zeros_like(lags)
stwoRANDarr=np.nan*np.zeros_like(lags)
for i in range(0,np.size(lags)):
print('Lag between', lags[i]-s_lag, 'and', lags[i]+s_lag)
good1, good2 = np.logical_and(dist >= lags[i]-s_lag, dist < lags[i]+s_lag).nonzero()
maskvec=mask.ravel()
Qvec=Qmap.ravel(); Qvec[maskvec > 1.]=np.nan
Uvec=Umap.ravel(); Uvec[maskvec > 1.]=np.nan
QvecRAND=QmapRAND.ravel(); QvecRAND[maskvec > 1.]=np.nan
UvecRAND=UmapRAND.ravel(); UvecRAND[maskvec > 1.]=np.nan
if np.logical_and(np.size(good1) > 0, np.size(good2) > 0):
# --------------------------------------------------------------------------------
Q1=Qvec[good1]
U1=Uvec[good1]
Q2=Qvec[good2]
U2=Uvec[good2]
good=np.logical_and(np.logical_and(np.isfinite(Q1),np.isfinite(Q2)),np.logical_and(np.isfinite(U1),np.isfinite(U2))).nonzero()
# From Eq. (11) in Planck XIX. A&A 576 (2015) A104
deltapsi=0.5*np.arctan2(Q1[good]*U2[good]-Q2[good]*U1[good], Q1[good]*Q2[good]+U1[good]*U2[good])
# From Eq. (6) in Planck XIX. A&A 576 (2015) A104
gooddeltapsi=deltapsi[np.isfinite(deltapsi).nonzero()]
weights=0.5*np.ones_like(gooddeltapsi)
stwo=np.sqrt(np.sum(weights*gooddeltapsi**2)/np.sum(weights))
# --------------------------------------------------------------------------------
Q1=QvecRAND[good1]
U1=UvecRAND[good1]
Q2=QvecRAND[good2]
U2=UvecRAND[good2]
good=np.logical_and(np.logical_and(np.isfinite(Q1),np.isfinite(Q2)),np.logical_and(np.isfinite(U1),np.isfinite(U2))).nonzero()
# From Eq. (11) in Planck XIX. A&A 576 (2015) A104
deltapsiRAND=0.5*np.arctan2(Q1[good]*U2[good]-Q2[good]*U1[good], Q1[good]*Q2[good]+U1[good]*U2[good])
# From Eq. (6) in Planck XIX. A&A 576 (2015) A104
gooddeltapsiRAND=deltapsiRAND[np.isfinite(deltapsiRAND).nonzero()]
weights=0.5*np.ones_like(gooddeltapsiRAND)
stwoRAND=np.sqrt(np.sum(weights*gooddeltapsiRAND**2)/np.sum(weights))
else:
print("No points in the selected lag range")
stwo=np.nan
stwoRAND=np.nan
stwoarr[i]=stwo
stwoRANDarr[i]=stwoRAND
#import pdb; pdb.set_trace()
return stwoarr, stwoRANDarr
# ------------------------------------------------------------------------------------------------------------------------------
def structurefunction(Qmap, Umap, lag=4.0, s_lag=1.0, mask=None, header=None):
""" Calculates the spatial correlation between im1 and im2 using the HOG method
Parameters
----------
Qmap : array corresponding to the first image to be compared
Umap : array corresponding to the second image to be compared
Returns
-------
hogcorr :
corrframe :
Examples
--------
"""
assert Qmap.shape == Umap.shape, "Dimensions of Qmap and Umap must match"
if mask is None:
mask=np.ones_like(Qmap)
assert Qmap.shape == mask.shape, "Dimensions of mask and Qmap must match"
Qmap[(mask==0.).nonzero()]=np.nan
Umap[(mask==0.).nonzero()]=np.nan
sz=np.shape(Qmap)
if header is None:
posx, posy = np.meshgrid(np.arange(0,sz[0]), np.arange(0,sz[1]))
else:
ra=header['CRVAL1']+header['CDELT1']*(np.arange(header['NAXIS1'])-header['CRPIX1'])
dec=header['CRVAL2']+header['CDELT2']*(np.arange(header['NAXIS2'])-header['CRPIX2'])
posx, posy=np.meshgrid(ra, dec)
x1, x2 =np.meshgrid(posx.ravel(), posx.ravel())
deltax=(x1-x2)
x1=None; x2=None
y1, y2 =np.meshgrid(posy.ravel(), posy.ravel())
deltay=(y1-y2)
y1=None; y2=None
print('Calculation of distances')
dist=np.sqrt(deltax**2+deltay**2)
deltax=None; deltay=None;
good1, good2 = np.logical_and(dist >= lag-s_lag, dist < lag+s_lag).nonzero()
stwo=np.nan
ngood=0.
maskvec=mask.ravel()
Qvec=Qmap.ravel(); Qvec[maskvec > 1.]=np.nan
Uvec=Umap.ravel(); Uvec[maskvec > 1.]=np.nan
# From Eq. (11) in Planck XIX. A&A 576 (2015) A104
if np.logical_and(np.size(good1) > 0, np.size(good2) > 0):
Q1=Qvec[good1]
U1=Uvec[good1]
Q2=Qvec[good2]
U2=Uvec[good2]
good=np.logical_and(np.logical_and(np.isfinite(Q1),np.isfinite(Q2)),np.logical_and(np.isfinite(U1),np.isfinite(U2))).nonzero()
deltapsi=0.5*np.arctan2(Q1[good]*U2[good]-Q2[good]*U1[good], Q1[good]*Q2[good]+U1[good]*U2[good])
ngood=np.size(good)
#deltapsi=0.5*np.arctan2(Qvec[good1]*Uvec[good2]-Qvec[good2]*Uvec[good1], Qvec[good1]*Qvec[good2]+Uvec[good1]*Uvec[good2])
#import pdb; pdb.set_trace()
# From Eq. (6) in Planck XIX. A&A 576 (2015) A104
gooddeltapsi=deltapsi[np.isfinite(deltapsi).nonzero()]
weights=0.5*np.ones_like(gooddeltapsi)
#stwo=np.std(deltapsi)
stwo=np.sqrt(np.sum(weights*gooddeltapsi**2)/np.sum(weights))
else:
print("No points in the selected lag range")
#print(stwo)
from scipy.stats import circstd
#print(circstd)
#import pdb; pdb.set_trace()
#return stwo
return {'S2': stwo, 'npairs': ngood}
# ------------------------------------------------------------------------------------------------------------------------------
def OLDstructurefunction(Qmap, Umap, lag=1.0, s_lag=0.5, mask=0.0, pitch=1.0):
sz=np.shape(Qmap)
x=np.arange(-sz[1]/2., sz[1]/2., 1.)
y=np.arange(-sz[0]/2., sz[0]/2., 1.)
xx, yy = np.meshgrid(x, y)
#index=np.arange(np.size(xx))
sfmap=0.*Qmap
valid=(mask > 0).nonzero()
for i in range(0, np.size(Qmap)):
if (mask[np.unravel_index(i,sz)] > 0):
diff=np.sqrt((xx[np.unravel_index(i,sz)]-xx)**2+(yy[np.unravel_index(i,sz)]-yy)**2)
good=np.logical_and(mask>0, np.logical_and(diff>lag-s_lag,diff<lag+s_lag)).nonzero()
#Qdiff=np.mean(Qmap[np.unravel_index(i,sz)]-Qmap[good])
#Udiff=np.mean(Umap[np.unravel_index(i,sz)]-Umap[good])
#sfmap[np.unravel_index(i,sz)]=0.5*np.arctan2(-Udiff, Qdiff)
Qdiff=Qmap[np.unravel_index(i,sz)]*Qmap[good]+Umap[good]*Umap[np.unravel_index(i,sz)]
Udiff=Qmap[np.unravel_index(i,sz)]*Umap[good]-Qmap[good]*Umap[np.unravel_index(i,sz)]
angles=0.5*np.arctan2(Udiff, Qdiff)
sfmap[np.unravel_index(i,sz)]=np.sum(angles**2)/float(np.size(angles)-1)
#print(sfmap[np.unravel_index(i,sz)])
#plt.imshow(sfmap*180./np.pi, origin='lower')
#plt.show()
goodsfmap=sfmap[valid]
deltapsi=anglemean(goodsfmap[np.isfinite(goodsfmap).nonzero()])*180./np.pi
print(deltapsi)
#import pdb; pdb.set_trace()
return deltapsi