New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add block reduce function for HDUs #88

Merged
merged 11 commits into from Mar 29, 2014

Conversation

Projects
None yet
3 participants
@ellisowen
Contributor

ellisowen commented Feb 26, 2014

Adds rebin_cubeHDU to gammapy/image/utils.py

A test case still needs to be added, but it is not clear to me how best to include this in the test_utils.py script.

@cdeil

This comment has been minimized.

Member

cdeil commented Feb 26, 2014

Try to fix the travis-ci test errors:

At the moment I see this error:
https://travis-ci.org/gammapy/gammapy/jobs/19658656#L720

ImportError: No module named skimage.measure

You need to mark the test as optional, because skimage is an optional dependency of astropy and gammapy. Here is information how to do it (or look around other test files in astropy and gammapy for examples).

@@ -830,3 +831,50 @@ def paste_cutout_into_image(total, cutout, method='sum'):
total.data[y : y + dy, x : x + dx] = cutout.data
else:
raise ValueError('Invalid method: {0}'.format(method))
def rebin_cubeHDU(image_HDU, factor, func=np.sum, cube=False):

This comment has been minimized.

@cdeil

cdeil Feb 26, 2014

Member

Python function names should use snake_case, i.e. in this case it should read rebin_cube_hdu.
And parameter names should also use snake case, i.e. it should be image_hdu instead of image_HDU.
(I know HDU is an abbreviation, but it should still be hdu here.)

@@ -6,7 +6,8 @@
from astropy.units import Quantity
from astropy.io import fits
from astropy.wcs import WCS
from skimage.measure import block_reduce

This comment has been minimized.

@cdeil

cdeil Feb 26, 2014

Member

You must move this import into the function, because skimage is an optional dependency.

#Put rebinned data into a fitsHDU
rebinned_image = fits.ImageHDU(data=image_max1, header=header)
return rebinned_image

This comment has been minimized.

@cdeil

cdeil Feb 26, 2014

Member

Add newline at end of file.

cdelt1 = header['CDELT1']
cdelt2 = header['CDELT2']
#Define new header values for new resolution
header['CDELT1'] = cdelt1*factor

This comment has been minimized.

@cdeil

cdeil Feb 26, 2014

Member

There should be spaces around operators, i.e. it should be cdelt1 * factor instead of cdelt1*factor.
To fix this, highlight all of the function in eclipse and then execute Source -> Format Code.

(You should get in the habit of doing this for all the code you write and commit, just don't mix formatting changes with real changes in a single commit.)

@cdeil

This comment has been minimized.

Member

cdeil commented Feb 26, 2014

Oh, the failed import problem was in gammapy/image/utils.py, not in gammapy/image/tests/test_utils.py as I thought.
But my comment still applies ... add a test and use the common syntax to mark this test as using skimage as an optional dependency.

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Feb 27, 2014

Updated - I think this addresses the issues you mentioned.

@cdeil

This comment has been minimized.

Member

cdeil commented Feb 27, 2014

I don't think the FITS header WCS parameters of your output HDU are currently correct.

Here's an example of similar code that computes new FITS header WCS parameters:
https://github.com/keflavich/FITS_tools/blob/master/FITS_tools/fits_overlap.py#L25

Start by taking one or two example images and running this fimgbin tool to rebin it and get a new FITS header WCS.

Then check how it modifies the header keywords by running ftlist and the Unit diff tool ... or use the ftdiff tool directly to print the difference in header keywords.

Introduces new test cases for calc_footprint and block_reduce_hdu, im…
…proves the code for block_reduce_hdu in gammapy/image/utils.py
@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Feb 28, 2014

I think the header keywords are OK. The check you suggested only leads to changes in NAXIS1, NAXIS2, CDELT1, CDELT2, CRPIX1, CRPIX2:
....
eowen@hfm-1307f:~/analyses/galactic_high_energy/data$ ftlist kathrin.fits
Print options: H C K I T [K] K
SIMPLE = T / conforms to FITS standard
BITPIX = -32 / number of bits per data pixel
NAXIS = 2 / number of data axes
NAXIS1 = 9400 / length of data axis 1
NAXIS2 = 500 / length of data axis 2
EXTEND = T
CDELT1 = -0.0199966430664062 / [Degrees]
CDELT2 = 0.0199999809265146 / [Degrees]
OBJECT = 'Unknown '
DATE = '2013-03-26'
GCOUNT = 1 / required keyword; must = 1
EQUINOX = 2000.0
XTENSION= 'IMAGE ' / IMAGE extension
CUNIT1 = 'deg '
CUNIT2 = 'deg '
CTYPE2 = 'GLAT-CAR'
ORIGIN = 'HESS Collaboration'
DATE-MAP= ' '
CTYPE1 = 'GLON-CAR'
REFERENC= 'Unknown '
EXTNAME = 'DIFFUSEMAX2D' / Extension name
BUNIT = 'None '
DATAMIN = 0.0
CRVAL2 = 0.0 / [Degrees]
CRPIX1 = 4700.0 / [Reference pixel: centre of the image]
CRPIX2 = 250.500000465662 / [Reference pixel: centre of the image]
CRVAL1 = 341.010000228882 / [Degrees]
TELESCOP= 'H.E.S.S.'
PCOUNT = 0 / required keyword; must = 0
DATAMAX = 2.64716035758283E-08
RADESYS = ' '
COMMENT [Type of sky map] = 'Unknown '
END
....

eowen@hfm-1307f:~/analyses/galactic_high_energy/data$ ftlist rebinned_kathrin2.fits
Print options: H C K I T [H] K
SIMPLE = T / file does conform to FITS standard
BITPIX = -32 / number of bits per data pixel
NAXIS = 2 / number of data axes
NAXIS1 = 4700 / length of data axis 1
NAXIS2 = 250 / length of data axis 2
EXTEND = T / FITS dataset may contain extensions
COMMENT FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H
CDELT1 = -3.99932861328124028555E-02 / [Degrees]
CDELT2 = 3.99999618530292033736E-02 / [Degrees]
OBJECT = 'Unknown '
DATE = '2013-03-26'
EQUINOX = 2000.0
CUNIT1 = 'deg '
CUNIT2 = 'deg '
CTYPE2 = 'GLAT-CAR'
ORIGIN = 'HESS Collaboration'
DATE-MAP= ' '
CTYPE1 = 'GLON-CAR'
REFERENC= 'Unknown '
BUNIT = 'None '
DATAMIN = 0.0
CRVAL2 = 0.0 / [Degrees]
CRPIX1 = 2.35025000000000000000E+03 / [Reference pixel: centre of the image]
CRPIX2 = 1.25500000232830998925E+02 / [Reference pixel: centre of the image]
CRVAL1 = 341.010000228882 / [Degrees]
TELESCOP= 'H.E.S.S.'
DATAMAX = 2.64716035758283E-08
RADESYS = ' '
COMMENT [Type of sky map] = 'Unknown '
HISTORY TASK: FIMGBIN on FILENAME: kathrin.fits
END
....

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Feb 28, 2014

Please can you review my last commit?

@coveralls

This comment has been minimized.

coveralls commented Feb 28, 2014

Coverage Status

Coverage remained the same when pulling 4aedd09 on ellisowen:image_rebinning into f87e7f6 on gammapy:master.

Array containing down-sampling integer factor along each axis.
func : function
Function object which is used to calculate the return value for each local block. This function must implement an axis parameter such as numpy.sum or numpy.mean.

This comment has been minimized.

@cdeil

cdeil Feb 28, 2014

Member

Break line into two or three (respect ~ 80 character line length limit).

Parameters
----------
image_hdu : `astropy.io.fits.ImageHDU`
Original Image HDU, unscaled

This comment has been minimized.

@cdeil

cdeil Feb 28, 2014

Member

No need for empty newline between listed parameters.

def calc_footprint(header):
"""Compute WCS corner positions.

This comment has been minimized.

@cdeil

cdeil Feb 28, 2014

Member

Please add Parameter and Returns section ... maybe even one example in the Examples section?

return rebinned_image
def calc_footprint(header):
"""Compute WCS corner positions.

This comment has been minimized.

@cdeil

cdeil Feb 28, 2014

Member

Add See also section linking to WCS.calcFootprint and mention that this one differs by using the real corners of the image.

rebinned_image = fits.ImageHDU(data=data_reduced, header=header)
return rebinned_image
def calc_footprint(header):

This comment has been minimized.

@cdeil

cdeil Feb 28, 2014

Member

Rename to wcs_footprint.

def block_reduce_hdu(input_hdu, factors, func=np.sum):
"""Merges pixels together to reduce the resolution of the image.
Sums contribution of merged pixels. Factor must be an integer.

This comment has been minimized.

@cdeil

cdeil Feb 28, 2014

Member

It's not really true that the function always sums ... it can do other things.
And its also not true that factor must be an integer ... it can be a list of integers.

Remove this description.
Reference the scikit-image function by putting it in single ticks.
Locally build the sphinx docs and check that the link in the html docs works.

@@ -830,3 +830,68 @@ def paste_cutout_into_image(total, cutout, method='sum'):
total.data[y : y + dy, x : x + dx] = cutout.data
else:
raise ValueError('Invalid method: {0}'.format(method))
def block_reduce_hdu(input_hdu, factors, func=np.sum):

This comment has been minimized.

@cdeil

cdeil Mar 1, 2014

Member

I think it would be better to make the API as similar to the scikit-image block_reduce as possible.
I.e. use the same parameter names (rename factors to block_size, add cval) ... you can even re-use parts of most parts of their docstring.

@cdeil

This comment has been minimized.

Member

cdeil commented Mar 1, 2014

This doesn't seem right ... I think the footprint should not change in this case, no?

In [1]: from astropy.io import fits

In [2]: hdu = fits.open('GalacticSurvey_Model_StdExt_02.fits')[0]

In [3]: print(hdu.shape) # note: ny, nx
(500, 9400)

In [4]: from gammapy.image import block_reduce_hdu, calc_footprint

In [5]: hdu_1_10 = block_reduce_hdu(hdu, (1, 10))

In [6]: print(hdu_1_10.shape)
(500, 940)

In [7]: calc_footprint(hdu.header)
Out[7]: 
array([[  74.9899979 ,   -4.98999989],
       [  74.9899979 ,    4.98999989],
       [ 247.0100021 ,    4.98999989],
       [ 247.0100021 ,   -4.98999989]])

In [8]: calc_footprint(hdu_1_10.header)
Out[8]: 
array([[ 74.9899979 ,  -4.80999989],
       [ 74.9899979 ,  94.98999788],
       [ 56.20999832,  94.98999788],
       [ 56.20999832,  -4.80999989]])

I think you have the order of the x and y axis wrong when computing the new header keywords.
This might help:
http://photutils.readthedocs.org/en/latest/photutils/index.html#coordinate-conventions-in-photutils

The axis order here should be defined by numpy, i.e. what scikit-image does with block_size.
There should be a test case that catches this!

@cdeil

This comment has been minimized.

Member

cdeil commented Mar 18, 2014

@ellisowen Could you have another look at this?

I'm still very keen on having this available, because it can speed up fitting and convolution a lot.

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Mar 18, 2014

Sure - I'll try and get this working this week (if not, will do it when I'm back in HD towards the end of next week).

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Mar 26, 2014

This seems to fix the issue - let me know if there is still a problem...

@@ -597,7 +597,7 @@ def solid_angle(image):
Returns
-------
area_image : `~astropy.units.quantity.Quantity`
Solid angle image (matching the input image) in steradians.
Solid angle image (matching the input imafrom astropy.wcs import WCSge) in steradians.

This comment has been minimized.

@cdeil

cdeil Mar 27, 2014

Member

I think you pasted from astropy.wcs import WCS here accidentally ... please remove.

@cdeil

This comment has been minimized.

Member

cdeil commented Mar 27, 2014

For the example above I still get a change in footprint when re-binning a survey map, which should not be the case:

In [3]: calc_footprint(hdu.header)
Out[3]: 
array([[  74.9899979 ,   -4.98999989],
       [  74.9899979 ,    4.98999989],
       [ 247.0100021 ,    4.98999989],
       [ 247.0100021 ,   -4.98999989]])

In [4]: calc_footprint(hdu_1_10.header)
Out[4]: 
array([[ 74.9899979 ,  -4.80999989],
       [ 74.9899979 ,  94.98999788],
       [ 56.20999832,  94.98999788],
       [ 56.20999832,  -4.80999989]])

Do you see this issue or does it work for you for a survey map?

Further debugging of the calc_footprint and block_reduce_hdu function…
…s in gammapy.image and updates test cases
@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Mar 27, 2014

This seems ok (there is a small change in the footprint for various rebinnings, but I don't think this can be avoided)

In [2]: image = fits.open('/home/eowen/analyses/galactic_high_energy/hess_maps/spec_hard_zeta_1_10TeV_STANDARD/final/derivedMaps_0.2.fits')[2]

In [3]: from gam
gammafit  gammapy   

In [3]: from gammapy.image import calc_footprint, block_reduce_hdu

In [4]: calc_footprint(image.header)
Out[4]: 
{'LOWER_LEFT': [array(247.01577949523966), array(-4.999995240941882)],
 'LOWER_RIGHT': [array(74.98422431945795), array(-4.999995240941882)],
 'TOP_LEFT': [array(74.98422431945795), array(4.99999522231542)],
 'TOP_RIGHT': [array(247.01577949523966), array(4.99999522231542)]}

In [5]: image_2 = block_reduce_hdu(image, (5, 10))

In [6]: calc_footprint(ima
image    image_2  

In [6]: calc_footprint(image_2.header)
Out[6]: 
{'LOWER_LEFT': [array(246.97578620910684), array(-4.909995326772565)],
 'LOWER_RIGHT': [array(74.94423103332514), array(-4.909995326772565)],
 'TOP_LEFT': [array(74.94423103332514), array(5.089995136484734)],
 'TOP_RIGHT': [array(246.97578620910684), array(5.089995136484734)]}
@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Mar 28, 2014

In [5]: image = make_empty_image(100, 100, proj='CAR')

In [6]: image_w = make_empty_image(101, 101, proj='CAR')

In [7]: image_1 = block_reduce_hdu(image, (10, 10))

In [8]: image_w_1 = block_reduce_hdu(image, (10, 10))
In [10]: calc_footprint(image.header)
Out[10]: 
{'LOWER_LEFT': [array(355.0), array(-5.0)],
 'LOWER_RIGHT': [array(5.0), array(-5.0)],
 'TOP_LEFT': [array(5.0), array(5.0)],
 'TOP_RIGHT': [array(355.0), array(5.0)]}

In [11]: calc_footprint(image_w.header)
Out[11]: 
{'LOWER_LEFT': [array(354.95), array(-5.050000000000001)],
 'LOWER_RIGHT': [array(5.050000000000001), array(-5.050000000000001)],
 'TOP_LEFT': [array(5.050000000000001), array(5.050000000000001)],
 'TOP_RIGHT': [array(354.95), array(5.050000000000001)]}

In [12]: calc_footprint(image_1.header)
Out[12]: 
{'LOWER_LEFT': [array(355.0), array(-5.0)],
 'LOWER_RIGHT': [array(5.0), array(-5.0)],
 'TOP_LEFT': [array(5.0), array(5.0)],
 'TOP_RIGHT': [array(355.0), array(5.0)]}
In [13]: calc_footprint(image_w_1.header)
Out[13]: 
{'LOWER_LEFT': [array(355.0), array(-5.0)],
 'LOWER_RIGHT': [array(5.0), array(-5.0)],
 'TOP_LEFT': [array(5.0), array(5.0)],
 'TOP_RIGHT': [array(355.0), array(5.0)]}

This seems to fix the problem (small correction to reference pixels was made to do this). Any non-integer pixel after the block reduce is applied is thrown away.

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Mar 28, 2014

Changes commited. Please review.

@cdeil

This comment has been minimized.

Member

cdeil commented Mar 28, 2014

Can you fix this travis-ci failure?
https://travis-ci.org/gammapy/gammapy/jobs/21779503#L1383

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Mar 29, 2014

Fixed!

cdeil added a commit that referenced this pull request Mar 29, 2014

@cdeil cdeil merged commit 9dacde7 into gammapy:master Mar 29, 2014

1 check passed

default The Travis CI build passed
Details
@cdeil

This comment has been minimized.

Member

cdeil commented Mar 29, 2014

Merged.
Thanks for your patience with this pull request!

I know most other code in gammapy doesn't have the same quality (good docs, tests),
but I hope this will improve over the next few months ...

@ellisowen ellisowen deleted the ellisowen:image_rebinning branch Jul 10, 2014

@cdeil cdeil changed the title from Image rebinning to Add block reduce function for HDUs Apr 8, 2015

@cdeil cdeil added the feature label Apr 8, 2015

@cdeil cdeil added this to the 0.1 milestone Apr 8, 2015

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment