Skip to content
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 methods for WCS maps #1281

Merged
merged 11 commits into from Feb 7, 2018

Conversation

Projects
None yet
3 participants
@woodmd
Copy link
Member

woodmd commented Feb 3, 2018

This PR implements some methods for WcsGeom and WcsNDMap that are inspired by methods of the same name in SkyImage:

  • crop/pad : Decrease/increase the size of a geometry/map by the given number of pixels.
  • downsample/upsample : Decrease/increase the resolution of a geometry/map by the given factor.

These methods operate only in the spatial dimension of the map. We could consider extending the methods to operate on non-spatial dimensions but this is probably a less common use case. Some parts of the current implementation are still missing -- namely support for all padding modes. pad and crop are also restricted to be symmetric (i.e. the same number of pixels added on either side of the image).

This PR also includes a new coadd method to Map which fills the contents of one map with another. For maps with the same WCS or HPX projection this method can be used to stack counts maps together even if they have different pixel sizes. Presently this method operates on the map in-place but we could also consider having it return a copy of the co-added map.

@cdeil cdeil self-assigned this Feb 4, 2018

@cdeil cdeil added the feature label Feb 4, 2018

@cdeil cdeil added this to the 0.8 milestone Feb 4, 2018

@cdeil
Copy link
Member

cdeil left a comment

@woodmd - Thanks!

I had a very quick look, left some inline comments.

@registerrier - Can you have a look, please?

@@ -265,6 +265,24 @@ def sum_over_axes(self):
"""Reduce to a 2D image by summing over non-spatial dimensions."""
pass

def coadd(self, map_in):
"""Fill this map with the contents of another map. This method can be

This comment has been minimized.

@cdeil

cdeil Feb 4, 2018

Member

I find "fill" here slightly confusing. Can you change the first line to more clearly say that this "adds" and doesn't "overwrite" the content?

This comment has been minimized.

@woodmd

woodmd Feb 4, 2018

Author Member

Elsewhere "fill" is used to mean increment the value (e.g. the fill_by methods). I can change to "Add the contents of map_in to this map."

This comment has been minimized.

@registerrier

registerrier Feb 4, 2018

Contributor

What is the use case of the coadd method for maps with different projections? Shouldn't it be possible only if geometries are aligned and bin size equal? Summing counts this way makes sense only if the final bin size is larger than that of map_in.

This comment has been minimized.

@woodmd

woodmd Feb 5, 2018

Author Member

The most common use-case would be for maps that have the same projection. In that case one can combine counts maps as long as the bin size of map_in is an integer multiple of the map bin size. If you're worried that people may misuse this method we could consider putting in some check to ensure the geometries are compatible.

npix = (self.npix[0] / factor, self.npix[1] / factor)
cdelt = (self._cdelt[0] * factor, self._cdelt[1] * factor)
wcs = get_resampled_wcs(self.wcs, factor, True)
return self.__class__(wcs, npix, cdelt=cdelt,

This comment has been minimized.

@cdeil

cdeil Feb 4, 2018

Member

I think it might be a good idea to add a copy method, and to always go through that copy method in other methods that make new, similar maps. The problem with the current pattern is that it's hard to maintain / be consistent when we add properties like unit or meta, it's easy to forget to update some of these self.__class__ callers.

A copy method that by default makes a copy of all data members in a single place is easier to maintain. So by default it copies everything, but if the caller passes in things they want differnt about the copy (like here wcs), then those are used to construct the new object.

I haven't fully thought this through, so I'm not sure it's a good suggestion that really helps. @woodmd - Thoughts?

This comment has been minimized.

@woodmd

woodmd Feb 4, 2018

Author Member

A copy method seems like it could be useful but I think requires some careful thought. Is this an idiom that you use elsewhere in gammapy or have you seen other packages use it? My feeling is that this would probably better fit in a separate PR where we can have a dedicated discussion on the topic.

This comment has been minimized.

@cdeil

cdeil Feb 4, 2018

Member

Is this an idiom that you use elsewhere in gammapy or have you seen other packages use it?

Not yet. It's something I wanted to introduce after seeing issues with methods that construct new objects in other places, e.g. in SkyImage / SkyCube.

My feeling is that this would probably better fit in a separate PR where we can have a dedicated discussion on the topic.

Agreed. You could leave a # TODO comment, although if you always use self.__class__ in gammapy.maps it's also easy to find all cases via text search.

@@ -355,17 +355,77 @@ def _reproject_hpx(self, geom, mode='interp', order=1):

return map_out

def pad(self, pad_width):
raise NotImplementedError
def pad(self, pad_width, mode='edge', cval=0):

This comment has been minimized.

@cdeil

cdeil Feb 4, 2018

Member

This pad method is pretty long and hard to read. I only had a very quick look and it isn't obvious to me where the default case mode='edge' is handled. maybe it would help to split this out into two methods that represent the if / else? I'm not sure it's useful to have the single pad that then dispatches to the two different algorithms, but if so, then the other two methods could be private (leading _). Also needs a docstring.

This comment has been minimized.

@woodmd

woodmd Feb 4, 2018

Author Member

Ok I can try to split this out into separate private methods. Regarding the docstring all of the abstract methods inherit their docstring from the Map class so that's why no docstring appears here.

@@ -875,6 +875,18 @@ def to_cube(self, axes):
return self.__class__(np.max(self.nside), self.nest, coordsys=self.coordsys,
region=self.region, conv=self.conv, axes=axes)

def pad(self, pad_width):
raise NotImplementedError

This comment has been minimized.

@cdeil

cdeil Feb 4, 2018

Member

Will there ever be pad & crop for the HEALPix classes? Remove these NotImplemented methods?

This comment has been minimized.

@woodmd

woodmd Feb 4, 2018

Author Member

Yes I'm working on them and could potentially add them to this PR. The logic is somewhat more complex since the "edge" of the map is not as well-defined.

return self.__class__(wcs, npix, cdelt=copy.deepcopy(self._cdelt),
axes=copy.deepcopy(self.axes))

def crop(self, crop_width):

This comment has been minimized.

@registerrier

registerrier Feb 4, 2018

Contributor

Why is the crop limited to a box around the center of the WcsGeom?

This comment has been minimized.

@woodmd

woodmd Feb 5, 2018

Author Member

Here the intent was to reproduce the functionality of the corresponding SkyImage methods. We could certainly have a discussion about extending the functionality. For what you're suggesting one approach might be a new method crop_to(geom) which would create a new map from the all pixels contained within the input geometry. Alternatively one could add an optional parameter to crop that allows you to shift the registration point of the box in X/Y.

@woodmd

This comment has been minimized.

Copy link
Member Author

woodmd commented Feb 5, 2018

I've now added crop and pad for HpxGeom and HpxNDMap. These methods remove/add N rings of pixels at the geometry edge. In the case of crop this means removing all pixels that have at least one neighbor not contained in the geometry while for pad it adds all pixels that have at least one neighbor inside the geometry. Since these algorithms work iteratively on rings they're obviously not efficient when dealing with a very large pad/crop width.

After experimenting with these there are some aspects of their behavior that's not completely ideal -- adding/removing a ring of pixels has the effect of distorting the shape of the geometry in unpredictable ways depending on where you are on the sky. If the pixel size is sufficiently small the distortion effect is pretty small but it's something to keep in mind. In light of this I was thinking that we might also want to an option to allow padding/cropping as a function of angle instead of pixel width. Thus the method signatures could look something like:

def crop(self, crop_width, unit='pix'):
    ...
def pad(self, pad_width, unit='pix'):
    ...

such that unit='pix' would preserve the current behavior while unit='deg' would pad with a number of pixels that encompass the chosen margin in angle. For circular HPX maps this would address the issue of preserving the shape when using pad/crop.

@cdeil cdeil merged commit 6099422 into gammapy:master Feb 7, 2018

2 checks passed

continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
@cdeil

This comment has been minimized.

Copy link
Member

cdeil commented Feb 7, 2018

Thanks!

@cdeil cdeil modified the milestones: 0.8, 0.7 Apr 11, 2018

@cdeil cdeil changed the title Implement methods for WCS maps Add methods for WCS maps Apr 11, 2018

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.