Converting between ICRS and Helioprojective observations

[Kamen Kozarev]

Hello,

I have some ground-based solar radio imaging observations, and would like to convert them to helioprojective coordinates, so I can compare them with space-based AIA/RHESSI/etc observations.

I am getting the following when I create a map from the input fits file:

print mwamap.coordinate_frame

<ICRS Frame>

/Users/kkozarev/anaconda2/lib/python2.7/site-packages/sunpy/map/mapbase.pyc:632: Warning: Missing metadata for heliographic latitude: assuming Earth-based observer
  [ 1.    ,  0.1875, -0.8125,  0.125 ,  0.3125],
/Users/kkozarev/anaconda2/lib/python2.7/site-packages/sunpy/map/mapbase.pyc:634: Warning: Missing metadata for Sun-spacecraft separation: assuming Sun-Earth distance
  [-0.6875, -0.3125,  0.8125,  0.0625,  0.1875],

Below is the header of the fits file. 
I'd appreciate your help.

Thanks,
Kamen







SIMPLE  =                    T /Standard FITS                                   
BITPIX  =                  -32 /Floating point (32 bit)                         
NAXIS   =                    4                                                  
NAXIS1  =                 1024                                                  
NAXIS2  =                 1024                                                  
NAXIS3  =                    1                                                  
NAXIS4  =                    1                                                  
EXTEND  =                    T                                                  
BSCALE  =   1.000000000000E+00 /PHYSICAL = PIXEL*BSCALE + BZERO                 
BZERO   =   0.000000000000E+00                                                  
BMAJ    =   7.185632493761E-02                                                  
BMIN    =   6.830347696940E-02                                                  
BPA     =   8.906390380859E+01                                                  
BTYPE   = 'Intensity'                                                           
OBJECT  = 'Sun     '                                                            
                                                                                
BUNIT   = 'Jy/beam '           /Brightness (pixel) unit                         
EQUINOX =   2.000000000000E+03                                                  
RADESYS = 'FK5     '                                                            
LONPOLE =   1.800000000000E+02                                                  
LATPOLE =  -1.518175136898E+01                                                  
PC01_01 =   1.000000000000E+00                                                  
PC02_01 =   0.000000000000E+00                                                  
PC03_01 =   0.000000000000E+00                                                  
PC04_01 =   0.000000000000E+00                                                  
PC01_02 =   0.000000000000E+00                                                  
PC02_02 =   1.000000000000E+00                                                  
PC03_02 =   0.000000000000E+00                                                  
PC04_02 =   0.000000000000E+00                                                  
PC01_03 =   0.000000000000E+00                                                  
PC02_03 =   0.000000000000E+00                                                  
PC03_03 =   1.000000000000E+00                                                  
PC04_03 =   0.000000000000E+00                                                  
PC01_04 =   0.000000000000E+00                                                  
PC02_04 =   0.000000000000E+00                                                  
PC03_04 =   0.000000000000E+00                                                  
PC04_04 =   1.000000000000E+00                                                  
CTYPE1  = 'RA---SIN'                                                            
CRVAL1  =   2.181370138385E+02                                                  
CDELT1  =  -5.555555555556E-03                                                  
CRPIX1  =   5.130000000000E+02                                                  
CUNIT1  = 'deg     '                                                            
CTYPE2  = 'DEC--SIN'                                                            
CRVAL2  =  -1.518175136898E+01                                                  
CDELT2  =   5.555555555556E-03                                                  
CRPIX2  =   5.130000000000E+02                                                  
CUNIT2  = 'deg     '                                                            
CTYPE3  = 'FREQ    '                                                            
CRVAL3  =   1.188400000115E+08                                                  
CDELT3  =   2.800000046235E+05                                                  
CRPIX3  =   1.000000000000E+00                                                  
CUNIT3  = 'Hz      '                                                            
CTYPE4  = 'STOKES  '                                                            
CRVAL4  =  -5.000000000000E+00                                                  
CDELT4  =   1.000000000000E+00                                                  
CRPIX4  =   1.000000000000E+00                                                  
CUNIT4  = '        '                                                            
PV2_1   =   0.000000000000E+00                                                  
PV2_2   =   0.000000000000E+00                                                  
RESTFRQ =   1.188367437383E+08 /Rest Frequency (Hz)                             
SPECSYS = 'TOPOCENT'           /Spectral reference frame                        
ALTRVAL =  -8.214682714364E+03 /Alternate frequency reference value             
ALTRPIX =   1.000000000000E+00 /Alternate frequency reference pixel             
VELREF  =                  259 /1 LSR, 2 HEL, 3 OBS, +256 Radio                 
COMMENT casacore non-standard usage: 4 LSD, 5 GEO, 6 SOU, 7 GAL                 
TELESCOP= 'MWA     '                                                            
OBSERVER= 'tfranzen'                                                            
DATE-OBS= '2015-11-04T03:43:45.250000'                                          
TIMESYS = 'UTC     '                                                            
OBSRA   =   2.181370138385E+02                                                  
OBSDEC  =  -1.518175136898E+01                                                  
OBSGEO-X=  -2.559454080476E+06                                                  
OBSGEO-Y=   5.095372146151E+06                                                  
OBSGEO-Z=  -2.849057186145E+06                                                  
DATE    = '2016-10-25T16:06:35.983000' /Date FITS file was written  

[Kamen Kozarev]

Anyone?

[Bin Chen]

Kamen,

I do not think there is an existing package in Sunpy that does the work of converting from celestial coordinates (RA and DEC) to helioprojective coordinates. You need to first find out the RA and DEC coordinate of the solar disk center as well as the p angle at your observation time. I usually use JPL Horizons, but python has a package, pyephem, to do this as well. Then, you figure out the offset of your map center in RA and DEC (do not forget the cos(dec) factor when converting d_RA to the actual arcsec values in projection), and rotate your image according to the p angle. If you wish to write out a fits file, you'd also need to update the tags that describe the coordinate system (CTYPE1 and CTYPE2 to HPLN_TAN and HPLT_TAN), pixel size/unit, etc.

The warning messages are about missing info of the observer. It does not hurt because MWA is certainly Earth based. To get rid of them, you'll need to add FITS tags HGLN_OBS and HGLT_OBS, and perhaps also RSUN_OBS and RSUN_REF. You may get the relevant information using sunpy.sun. 

Regards,

Bin

[Jack Ireland]

Kamen, Bin,

There is code in SunPy 0.8 that I think will do what you want - superimpose ground based images with space based.  It was developed for composing eclipse photos with space based data.

Please see

https://www.numfocus.org/blog/how-to-compare-photos-of-the-solar-eclipse-using-python-sunpy/

https://github.com/sunpy/solar-eclipse/blob/master/notebooks/Processing%20your%20Eclipse%20Photo%20with%20SunPy.ipynb

https://github.com/sunpy/solar-eclipse

Thanks,

Jack

.

--
You received this message because you are subscribed to the Google Groups "SunPy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sunpy+unsubscribe@googlegroups.com.
To post to this group, send email to su...@googlegroups.com.
Visit this group at https://groups.google.com/group/sunpy.
For more options, visit https://groups.google.com/d/optout.

[Jack Ireland]

Also useful is

http://docs.sunpy.org/en/stable/code_ref/coordinates.html#transforming-between-coordinate-frames

Jack

[Kamen Kozarev]

Thanks, Bin and Jack! 

To unsubscribe from this group and stop receiving emails from it, send an email to sunpy+un...@googlegroups.com.

[Bin Chen]

Thanks, Jack! The links look very useful. Has anyone tried to use Skycoord for a direct transformation from celestial ICRS to helioprojective? 

Kamen, I'd appreciate it if you could let me know if that works!

-Bin

To unsubscribe from this group and stop receiving emails from it, send an email to sunpy+un...@googlegroups.com.

[Jack Ireland]

Bin,

The section of the notebook

https://github.com/sunpy/solar-eclipse/blob/master/notebooks/Processing%20your%20Eclipse%20Photo%20with%20SunPy.ipynb

about overplotting Regulus on an eclipse photo has an example of ICRS to helioprojective.

Jack

To unsubscribe from this group and stop receiving emails from it, send an email to sunpy+unsubscribe@googlegroups.com.

[Bin Chen]

Got it. I'll give it try. Thanks much!

[Shih, Albert Y. (GSFC-6710)]

Hi, Bin and others,

     I wrote the code that enables the transformations from Astropy’s frames to SunPy’s frames, so let me know if there are any issues.  I’ll point out some potential sticking points:

 

·         You must provide distance information in addition to RA/dec to convert from ICRS to helioprojective due to the shift in origin.  (The distance can be a faked large value, of course.)

·         Transforming in the reverse direction – helioprojective to most of Astropy’s frames – currently only works for coordinates that point to something on the solar surface.

Albert

[Bin Chen]

Thanks, Albert. I tried the codes that Jack pointed out and ploted RA and DEC of the solar disk center on to an AIA map. But there appears to be some offset of ~30" in solar Y if I assume the distance is very large (1 ly). The offset would be much larger if I assume the distance is 1 AU. Here is a notebook that describes my procedure and results:

https://github.com/binchensun/suncasa/blob/master/notebooks/ICRS_on_AIA.ipynb. Any clue?

This is nice for over plotting something in ICRS coordinates to a map with helioprojective coordinates. But I guess what Kamen (and other solar radio astronomers) wants is to convert a map originally in ICRS to helioprojective and update the image/header directly. I wrote something in a more or less manual way (as I described in my first message) to convert a CASA image in ICRS to a FITS file with helioprojective coordinates: https://github.com/binchensun/suncasa/blob/master/utils/helioimage2fits.py (caution: it requires CASA and may not be very reader friendly!). But there could be a better way of doing this by using WCS, I suppose.

Regards,

Bin

[Shih, Albert Y. (GSFC-6710)]

Hi, Bin,

     Ah, you made a mistake when using the output of JPL Horizons.  Since you obtained the ICRF RA/dec of the Sun for a geocentric observer, you can’t use ICRS, because the origin is not the solar-system barycenter.  So, instead of creating the SkyCoord using the ICRS frame, the (more) correct frame to use is the GCRS frame.

     However, I still get a ~15-arcsec mismatch, primary in helioprojective X.  I believe this is because the GCRS frame factors in aberration (due to observer motion), while ICRF coordinates do not.  I’m not immediately sure what the best way is to insert ICRF coordinates with a non-barycenter origin.

     For this particular test case, of course, you can just use Astropy/SunPy to get the position of the Sun.  At this time, the Sun in GCRS (which, again, includes aberration) is at:

-----

>>> sun = SkyCoord(0*u.deg, 0*u.deg, 0*u.AU, obstime='2014-11-01 19:15:30', frame='heliographic_stonyhurst')

>>> sun.gcrs.to_string('hmsdms')

'14h26m45.5162s -14d31m16.1748s'

-----

[Bin Chen]

Thank you, Albert!

-Bin

[Kamen Kozarev]

Hello everyone,

I think I am getting close to generating a SunPy Map for MWA radio data. However, it seems that I need to 'flip' the vertical axis of the MWA map - Is there a way to do this with sunpy? Also, if anyone has general comments on this, I'd be grateful. I've used the solar eclipse code as well.

The Python notebook is here:

https://github.com/kkozarev/notebooks/blob/master/test_sunpy_mwa_demo.ipynb

Kamen

On Tuesday, November 7, 2017 at 2:24:11 PM UTC+2, Kamen Kozarev wrote:

[Kamen Kozarev]

Dear Albert,

I think I am getting close to generating a SunPy Map for MWA solar radio observations. However, it seems that I need to 'flip' the vertical axis of the MWA map - Is there a way to do this with sunpy? Also, if you have general comments on this, I'd be grateful. I've used the solar eclipse code as well.

The Python notebook is here:

https://github.com/kkozarev/notebooks/blob/master/test_sunpy_mwa_demo.ipynb

Kamen

[Jack Ireland]

Kamen,

If I understand this correctly, the problem is that the composite map flips the MWA map upside down?

Jack

To unsubscribe from this group and stop receiving emails from it, send an email to sunpy+unsubscribe@googlegroups.com.

[Shih, Albert Y. (GSFC-6710)]

Hi, Kamen,

     I’m not familiar with MWA files, but I have some comments towards what I think is going on:

 

·         The MWA file appears to have CDELT values indicating a left-handed coordinate system, which is presumably why the plot of the MWA map is using an X axis that goes from positive at the left to negative at the right.

·         You appear to be attempting to de-rotate the MWA map, but its PC matrix appears to be the identity (i.e., no rotation).

·         You appear to be attempting to de-rotate the MWA map by the AIA map’s rotation, which looks like it’s also the identity.

·         Your de-rotated map – which, as per above, should in theory be unchanged – may have mixed weirdly with its left-handed coordinate system, which is why it’s now the Y axis that is “reversed” rather than the X axis.

·         Given that all of the plots appear to be consistent, just with one axis or the other reversed, I hope that the composite map is showing the MWA contours in the correct place relative to the AIA map.

 

I recommend avoiding the left-handed coordinate system at the outset.  When you read in the MWA file, I’d flip the data array – numpy.fliplr()? – and removing the negative sign on the first element of CDELT.

 

Albert

[Bin Chen]

You are right, Albert. The X axis in MWA maps is most likely right ascension, which increases toward the east, opposite from HPLN-TAN (solar X).

-Bin

To unsubscribe from this group and stop receiving emails from it, send an email to sunpy+unsubscribe@googlegroups.com.

To post to this group, send email to su...@googlegroups.com.
Visit this group at https://groups.google.com/group/sunpy.
For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to the Google Groups "SunPy" group.

To unsubscribe from this group and stop receiving emails from it, send an email to sunpy+unsubscribe@googlegroups.com.

To post to this group, send email to su...@googlegroups.com.
Visit this group at https://groups.google.com/group/sunpy.
For more options, visit https://groups.google.com/d/optout.

--

Dr. Bin Chen
Assistant Professor of Physics

Center for Solar-Terrestrial Research

New Jersey Institute of Technology

101 Tiernan Hall

323 Dr. Martin Luther King Blvd, Newark, NJ 07102

(973) 596-3565

bin....@njit.edu

https://web.njit.edu/~binchen

[Shih, Albert Y. (GSFC-6710)]

Hi, Bin,

     Ah, quite right; I missed that the MWA header does in fact say that the coordinate system is RA/dec.

 

Kamen,

     Given that, it looks like you’ve incorrectly assigned the helioprojective coordinate system to the MWA map, so ignore the recommendation at the end of my email.  You need to fix the coordinate-system issue.  Presumably the MWA contours are showing up in the wrong place relative to the AIA map.

 

Albert

-Bin

Dr. Bin Chen
Assistant Professor of Physics

Center for Solar-Terrestrial Research

New Jersey Institute of Technology

 

101 Tiernan Hall

323 Dr. Martin Luther King Blvd, Newark, NJ 07102

(973) 596-3565

bin....@njit.edu

https://web.njit.edu/~binchen

--

[Kamen Kozarev]

Well, I tried to follow the example in the solar eclipse notebook.

If I set the wcs header as follows:

fh = pyfits.open(mwafile)

data=fh[0].data[0,0,:,:]

header = fh[0].header

solar_rotation_angle = sunpy.coordinates.get_sun_orientation(loc, header['date-obs'])

w = wcs.WCS(naxis=2)

w.wcs.dateobs = header["DATE-OBS"]

w.wcs.crpix = [header['crpix1'],header['crpix2']]

w.wcs.ctype = [header['ctype1'],header['ctype2']]

w.wcs.cunit = ['deg', 'deg']

w.wcs.crval = [header['crval1'],header['crval2']]

w.wcs.cdelt = [header['cdelt1'],header['cdelt2']]

#Convert to a header which will be used to generate the map

mapheader = dict(w.to_header())

mapheader.update({'CROTA': solar_rotation_angle.to('deg').value})

mapheader.update({'DSUN_OBS': dsun.to('m').value})

mapheader.update({'RSUN': dsun.to('m').value})

mapheader.update({'RSUN_OBS': np.arctan(sunpy.sun.constants.radius / dsun).to('arcsec').value})

mwamap = Map(data, mapheader)

,I get an error when trying to convert to helioprojective:

aiamap = Map(aiafile)

fig = plt.figure(figsize=(10,10))

ax = plt.subplot(projection=mwamap)

mwamap.plot(axes=ax)

mwamap.draw_grid(axes=ax)

mwamap.draw_limb(axes=ax)

plt.show()

UnitsError: The input ICRS coordinates do not have length units. This probably means you created coordinates with lat/lon but no distance.  Heliocentric<->ICRS transforms cannot function in this case because there is an origin shift.

Any thoughts?