3

I have a set of data given here where in the first and second columns there are the sky coordinates (ra,dec), respectively and in the third and forth, the coordinates in a Cartesian system (x,y).

enter image description here

I need to make a two-dimensional interpolation surface using coordinates x and y and another using Ra and Dec. The problem is the existence of masked regions, as shown in the figure above. I can illustrate the missing data just by plotting them (There is non NaN value in the catalogue). That is what I so far tried and didn't give the right answer:

from scipy.interpolate import griddata
import numpy as np
import matplotlib.pyplot as plt

data = np.loadtxt('test.asc')

ra = data[:,0]
dec = data[:,1]
Xpos = data[:,2]
Ypos = data[:,3]

xi = np.linspace(Xpos.min(), Xpos.max(), 1000)
yi = np.linspace(Ypos.min(), Ypos.max(), 1000)
xi, yi = np.meshgrid(xi, yi, copy=False)
ra_int = griddata(data[:,2:4], ra, (xi.flatten(), yi.flatten()),
                  method='cubic')
dec_int = griddata(data[:,2:4], dec, (xi.flatten(), yi.flatten()),
                   method='cubic')

Using griddata fails and return just NaN values. Is there any way to do this interpolation in order to estimate the values of Ra and Dec from a given x and y coordinates even in the masked regions (map from x and y to ra and dec)?

Dalek
  • 4,168
  • 11
  • 48
  • 100
  • You might try to fill the voids as shown in this post http://stackoverflow.com/questions/20753288/filling-gaps-on-an-image-using-numpy-and-scipy/20902981#20902981 and interpolate thereafter. – Christian K. Sep 02 '14 at 13:32
  • Actually I did not want to point you to that specific answer (and I cannot edit the comment), but the problem is similar to yours. I guess Alex I's answer can suite very well. – Christian K. Sep 02 '14 at 13:39
  • @ChristianK. It seems very intriguing but then how could I map from one plane to the other plane? It is not very clear for me how it should be done! – Dalek Sep 02 '14 at 13:42
  • It seems I misunderstood your question. I was assuming that part of the problem was interpolation of a z = f(x,y) surface which has gaps, but this is obviously not the case. Actually I have no idea what you intend to do. `griddata` is for interpolation of f(x1,x2,x3,...,xn) type data. You seem to want to do soemthing else... – Christian K. Sep 02 '14 at 17:03
  • The RA and Dec are in degrees I assume (equatoreal coordinate system commonly used in astronomy data) but what exactly are the cartesian set of coordinates ? (some projection onto Earth surface at specific geolocation ?) without this info is hard to answer because we do not know the dependency between these two spaces (most likely non linear and not just mirrored and scaled ...) – Spektre Sep 03 '14 at 06:58
  • @Spektre The `X` and `Y` are the pixel coordinate of objects on the CCD image. – Dalek Sep 03 '14 at 07:41
  • in that case it is non linear (but not too much curved) ... just solve the equation for projection from spherical surface onto plane and calibrate the constants from midle and edge points in your dataset – Spektre Sep 03 '14 at 09:01
  • @Spektre can you elaborate a bit more what do you mean by *calibrate the constants from midle and edge points in your dataset*? – Dalek Sep 03 '14 at 09:12
  • added more info as answer it should work if the CCD is not axis aligned to equator then you need to apply also rotation by Z axis on the final result – Spektre Sep 03 '14 at 09:24
  • edge points are any near corners of the CCD and middle point is any point near middle of CCD – Spektre Sep 03 '14 at 09:27

1 Answers1

1

If I get it right then it is like this:

CCD projection

just shift the Cartesian coordinate system to middle of the CCD and also the Equatoreal coordinates to middle of CCD. Then compute x,y separately. The only thing you need is to compute focus length f separately for x and y !!!


pos is the cartesian coordinate (x or y)
ang is the equatoreal coordinate (RA or Dec)

  1. get edge point from the database

    shift the angles to middle of CCD

  2. compute focus (fx,fy) from it

    f = pos/tan(ang)

  3. now you can compute the projection for any entry in dataset

    shift the angles to middle of CCD then compute x,y by

    pos=f*tan(ang)

    shift back from CCD middle to original Cartesian coordinates. You should check few points if is this approach correct

[notes]

x axis is mirrored in your output so just use x=-x at the end before shifting back to original Cartesian coordinates or leave focus f negative.

if your CCD is not axis aligned to equator then you need to compute the rotation (angle between X axis and equator) and apply rotation around Z axis after conversion before shifting back...

Spektre
  • 49,595
  • 11
  • 110
  • 380