-
Notifications
You must be signed in to change notification settings - Fork 0
/
simple_map_xy.pro
148 lines (121 loc) · 3.05 KB
/
simple_map_xy.pro
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
PRO simple_map_xy, data, x, y, map, nsamp, weight, weight_map, flags, flags_map
;+
; NAME:
; SIMPLE_MAP
;
;
; PURPOSE:
; Takes a series of data samples with coordinates and maps them to
; an image header
;
;
; CATEGORY:
;
;
;
; CALLING SEQUENCE:
; SIMPLE_MAP, data, x, y, map, [weight, weight_map]
;
;
; INPUTS:
; data: data samples
; x: coordinates of samples
; y: coordinates of samples
;
;
;
; OPTIONAL INPUTS:
; weight: weights for data sample
; flags: flags (or other variable you want co-added)
;
; KEYWORD PARAMETERS:
;
;
;
; OUTPUTS:
; map: map with geometry defined by header and data-samples co-added
; at nearest pixel
; nsamp: number of sampels contributing to a map pixel (output)
;
;
; OPTIONAL OUTPUTS:
; weight_map: weight map on same system as map
; flags_map: flags map co-added on same system as map
;
;
; COMMON BLOCKS:
;
;
;
; SIDE EFFECTS:
;
;
;
; RESTRICTIONS:
;
;
;
; PROCEDURE:
; uses the histogram function for speed.
;
;
; EXAMPLE:
;
;
;
; MODIFICATION HISTORY:
;
; Original version December 19th 2007. Seb Oliver
;
;-
; extracting dimension of image from header
IF NOT keyword_set(naxis1) THEN BEGIN
IF keyword_set(map) THEN BEGIN
sz = size(map)
naxis1 = sz[1]
naxis2 = sz[2]
ENDIF ELSE BEGIN
message, 'naxis1 not defined'
ENDELSE
ENDIF
IF NOT keyword_set(naxis2) THEN naxis2 = naxis1
map = fltarr(naxis1, naxis2)
weight_map = fltarr(naxis1, naxis2)
flag_map = fltarr(naxis1, naxis2)
IF n_params() LE 6 THEN BEGIN
weight = replicate(1, n_elements(data))
DO_weight = 0
ENDIF ELSE DO_weight = 1
IF n_params() LE 8 THEN DO_flag = 0 ELSE DO_flag = 1
; weighting data samples
data_tmp = data*weight
; checking which samples within field
inmap = where(x GT -0.5 AND x LT naxis1-0.5 AND y GT -0.5 AND y LT naxis2-0.5, ngood)
IF ngood LE 0 then message, 'No good data'
; rounding to nearest pixel
ix = round(x[inmap])
iy = round(y[inmap])
; creating 1D index
index = long(iy)*naxis1+long(ix)
; using histogram to map samples to image space
nsamp = histogram(index, reverse_indi=r, min=0, max=long(naxis1)*naxis2-1)
; only work with map pixels that have some data
good = where(nsamp GT 0, ngood)
IF ngood LE 0 THEN message, 'No good data'
; this loop appears to be necessary and is the shortest loop I can
; imagine, being only 1
FOR j=0L, n_elements(good)-1L DO BEGIN
i = good[j]
; IF nsamp[i] EQ 10 THEN stop
; See IDL manual to see where next two lines comes from
IF R[i] NE R[i+1] THEN MAP[i] = total(data_tmp[inmap[ [R[R[I] : R[i+1]-1]]]] ) $
ELSE print, 'This should not happen'
IF keyword_set(DO_weight) THEN weight_MAP[i] = total(weight[inmap [R[R[I] : R[i+1]-1]]] )
IF keyword_set(DO_flags) THEN flag_MAP[i] = total(flags[inmap [R[R[I] : R[i+1]-1]] ])
ENDFOR
;reformating arrays to 2D
map = reform(map, naxis1, naxis2)
nsamp = reform(nsamp, naxis1, naxis2)
IF keyword_set(DO_weight) THEN weight_map = reform(weight_MAP, naxis1, naxis2)
IF keyword_set(DO_flags) THEN weight_map = reform(flags_MAP, naxis1, naxis2)
END