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 rotation support to the HEALPix projection #1638

Merged
merged 10 commits into from Oct 1, 2019

Conversation

Schneegans
Copy link
Contributor

@Schneegans Schneegans commented Sep 25, 2019

Motivation

At the German Aerospace Center (DLR) we are developing a VR-Framework called CosmoScout VR. This software can render highly detailed level-of-detail planetary bodies in real-time. The required surface tiles are downloaded from a MapServer instance and are provided in the HEALPix projection. The advantage is, that the data density is equal all around the globe without any singularities at the poles.

However, in order to retrieve the HEALPix base patches as image tiles, the map has to be rotated by 45°. Here is an image from our documentation which illustrates this with the figures in the center row:

image

Implemented Changes

The file src/projections/healpix.cpp has been changed to accept an optional +rot= parameter, which has to be given in degrees and defaults to zero.

Usage

We use this in a custom projection like this and it works very well:

<900914> +proj=healpix +lon_0=0 +x_0=2.5 +y_0=2.5 +a=0.900316316157106 +rot=45 +no_defs <>

Copy link
Member

@kbevers kbevers left a comment

Very interesting. I can't see why we shouldn't include this. Before we do that though you need to update the documentation. Start here https://github.com/OSGeo/PROJ/blob/master/docs/source/operations/projections/healpix.rst.

Also, does it make sense to add the rotation parameter to the rHEALPix projection as well?

@@ -78,6 +79,12 @@ static double sign (double v) {
return v > 0 ? 1 : (v < 0 ? -1 : 0);
}

static PJ_XY rotate(PJ_XY p, double a) {
Copy link
Member

@kbevers kbevers Sep 25, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I suggest you use angle instead of a here. Makes it immediately clear what's going on

@@ -622,6 +637,8 @@ PJ *PROJECTION(healpix) {
P->opaque = Q;
P->destructor = destructor;

Q->rot = pj_param(P->ctx, P->params,"drot").f * M_PI / 180.0;
Copy link
Member

@kbevers kbevers Sep 25, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Improve readability a bit:

Suggested change
Q->rot = pj_param(P->ctx, P->params,"drot").f * M_PI / 180.0;
double angle = pj_param(P->ctx, P->params,"drot").f;
Q->rot = PJ_TORAD(angle);

@rouault
Copy link
Member

@rouault rouault commented Sep 25, 2019

One thought I add is that such a +rot parameter could be something generic to all projections, such as +x_0, +y_0. But we could potentially move that to generic code at a later stage if we find that several projections need this.
Was wondering if rot_xy might not be more expressive about in which coordinate space the rotation is done ?

@rouault
Copy link
Member

@rouault rouault commented Sep 25, 2019

Extending heapix tests in test/gie/builtins.gie to test this new parameter would also be welcome

@kbevers
Copy link
Member

@kbevers kbevers commented Sep 25, 2019

One thought I add is that such a +rot parameter could be something generic to all projections, such as +x_0, +y_0. But we could potentially move that to generic code at a later stage if we find that several projections need this.

Like a rotation operation? Something like

+proj=pipeline +step +proj=healpix ... +step +proj=rot +angle=45.0

perhaps?

@Schneegans
Copy link
Contributor Author

@Schneegans Schneegans commented Sep 25, 2019

Thanks for the feedback! I'll update the PR as soon as possible.

Also, does it make sense to add the rotation parameter to the rHEALPix projection as well?

I like the idea of supporting this parameter for all projections, however for the rHEALPix projection I don't see any particular use case.

Extending heapix tests in test/gie/builtins.gie to test this new parameter would also be welcome

I will do this. Are these tests auto-generated? Or are the numbers somehow computed by hand? Or shall I just use the current implementation and assume that the values are correct?

@kbevers
Copy link
Member

@kbevers kbevers commented Sep 25, 2019

Or shall I just use the current implementation and assume that the values are correct?

Presumably you have tested elsewhere that your implementation actually works. If that is the case, then yes, just use values generated by the new code in the tests. They mainly exist as regression tests. When test material from a difference source is available, like test coordinates from a scientific paper, we like to use those as well but that is not possible in this case.

@rouault
Copy link
Member

@rouault rouault commented Sep 25, 2019

Like a rotation operation? Something like

yeah, I was thinking to that. Could be useful as an operation, but the issue is as soon as you define a pipeline, this is no longer a CRS definition, and you'll struggle using that in a Mapserver context where you define a CRS for your project or layer

Are these tests auto-generated? Or are the numbers somehow computed by hand? Or shall I just use the current implementation and assume that the values are correct?

manually. Just add a new entry and run "make check". It will give you the output result it found if not consistent with the one you put, so you can bootstrap with dummy output values, and use the ones it will get to you. But you're supposed to do some validation of your input/output tuples ;-)

@kbevers
Copy link
Member

@kbevers kbevers commented Sep 25, 2019

yeah, I was thinking to that. Could be useful as an operation, but the issue is as soon as you define a pipeline, this is no longer a CRS definition, and you'll struggle using that in a Mapserver context where you define a CRS for your project or layer

I see that as a Mapserver problem :-) But seriously, this could be dealt with by setting up a WKT string that includes a rotation step, no? Provided that the rotation operation is put in the WKT machinery of course.

@rouault
Copy link
Member

@rouault rouault commented Sep 25, 2019

I see that as a Mapserver problem :-)

yeah, which means my problem at some point (I've a port of MapServer to PROJ 6 API almost ready, and this wasn't an easy ride. And MapServer will still continue to use heavily PROJ strings for CRS in a foreseable future)

this could be dealt with by setting up a WKT string that includes a rotation step, no?

Yes, likely a DerivedProjectedCRS with a deriving conversion being a rotation. They will have derived projected CRS in EPSG v10 (see https://www.iogp.org/blog/epsg/upgrade-of-epsg-dataset-data-model/), but I haven't looked if the EPSG guidance note 7-2 describes a rotation method. Probably. Actually the +proj=affine ( https://proj.org/operations/transformations/affine.html ) could already be used for that (not completely sure to remember why it's classified as a transformation rather than a conversion, but whatever)


To come back to this pull request, what is the rotation center: (0, 0) or (x_0, y_0) ?

@busstoptaktik
Copy link
Member

@busstoptaktik busstoptaktik commented Sep 25, 2019

@kbevers said:

Like a rotation operation? Something like

We have that already, I believe, as a special case of the 4 parameter Helmert

@kbevers
Copy link
Member

@kbevers kbevers commented Sep 25, 2019

We have that already, I believe, as a special case of the 4 parameter Helmert

Correct. Although you have to give the rotation in arc seconds

@Schneegans
Copy link
Contributor Author

@Schneegans Schneegans commented Sep 25, 2019

To come back to this pull request, what is the rotation center: (0, 0) or (x_0, y_0) ?

I think the correct answer to this question is (x_0, y_0).

The first figure in the second row of the picture I included above shows a result of +x_0=2.5 +y_0=2.5 +rot_xy=45 (plus some scaling).

test/gie/builtins.gie Show resolved Hide resolved
@Schneegans
Copy link
Contributor Author

@Schneegans Schneegans commented Oct 1, 2019

Are there any updates for this one? I hope that I resolved all suggested changes 😄

@rouault
Copy link
Member

@rouault rouault commented Oct 1, 2019

Looks good to me. I let @kbevers do his check. (we should use "squash and merge" to fold the commits)

@kbevers kbevers merged commit 50a1821 into OSGeo:master Oct 1, 2019
4 checks passed
@kbevers
Copy link
Member

@kbevers kbevers commented Oct 1, 2019

Looks good to me

@kbevers kbevers added this to the 7.0.0 milestone Oct 1, 2019
@Schneegans Schneegans deleted the feature/healpix-rotation branch Oct 1, 2019
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

4 participants