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

added CRS mapping to/from CF 1.8; added to_proj4_dict #244

Merged
merged 3 commits into from
Apr 11, 2019

Conversation

snowman2
Copy link
Member

@snowman2 snowman2 commented Apr 6, 2019

Addresses #226

This will need a good review.

I am unsure about what the CF unit mapped to. It is either units or vunits. I am guessing vunits at the moment, but it is only a guess.

Others I think are just another maping to datum, but not sure. Maybe there is a vdatum in PROJ I am missing.

  • geoid_name (OGC WKT VERT_DATUM - ex GEOID12B)
  • geopotential_datum_name (OGC WKT VERT_DATUM - ex NAVD88)

@snowman2
Copy link
Member Author

snowman2 commented Apr 6, 2019

datum is just horizontal I believe.

@snowman2
Copy link
Member Author

snowman2 commented Apr 6, 2019

I am learning towards unit mapping to units as it seems the z coordinate will also have units, and that would correspond to the vunits I would think.

@djhoese
Copy link
Contributor

djhoese commented Apr 7, 2019

@snowman2 This is amazing. I think some good people to look at this would be @dopplershift, @mraspaud, and @pnuu. I'll try to find time to review this in the next couple days.

@jswhit
Copy link
Collaborator

jswhit commented Apr 8, 2019

Looks good. I just wonder if a warning should be issued if a parameter doesn't map to CF (instead of silently ignoring them).

Copy link
Contributor

@djhoese djhoese left a comment

Choose a reason for hiding this comment

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

I had a couple questions that I commented on specifically. Otherwise, my other question is whether it would be better to not copy the in_cf dictionary in from_cf but rather let KeyErrors handle the pop'd keys in the final for loop.

Again, this is great.

pyproj/crs.py Outdated Show resolved Hide resolved
pyproj/crs.py Outdated Show resolved Hide resolved
@snowman2
Copy link
Member Author

snowman2 commented Apr 8, 2019

The popping was from another iteration where I was rasing an exception when a parameter didn't match. I think @jswhit has a valid point about warning the user. If we do warn the user, the current copy/pop implementation should stay. If not, it can be changed.

I added warnings in the docs, but those aren't always read. So, putting it in the code will make it more likely that the user will notice.

@snowman2
Copy link
Member Author

snowman2 commented Apr 9, 2019

It seems there is a warning raised almost every single time the function to_cf is called.

  /home/snowal/scripts/pyproj/pyproj/crs.py:362: UserWarning: PROJ parameters not mapped to CF: ('k', 'no_defs', 'type')
    "PROJ parameters not mapped to CF: {}".format(tuple(skipped_params))

t /home/snowal/scripts/pyproj/pyproj/crs.py:362: UserWarning: PROJ parameters not mapped to CF: ('no_defs', 'type')
    "PROJ parameters not mapped to CF: {}".format(tuple(skipped_params))

test/test_crs_cf.py::test_cf_from_utm
  /home/snowal/scripts/pyproj/pyproj/crs.py:362: UserWarning: PROJ parameters not mapped to CF: ('zone', 'no_defs', 'type')
    "PROJ parameters not mapped to CF: {}".format(tuple(skipped_params))

test/test_crs_cf.py::test_cf_rotated_latlon
  /home/snowal/scripts/pyproj/pyproj/crs.py:362: UserWarning: PROJ parameters not mapped to CF: ('o_proj', 'type')
    "PROJ parameters not mapped to CF: {}".format(tuple(skipped_params))

That is kind of annoying, but good to know as well. So, I am not sure how I feel about it.

Thoughts?

@djhoese
Copy link
Contributor

djhoese commented Apr 9, 2019

My initial reaction is to mark certain parameters as "can't be converted" by storing them with None values in the dictionary and detecting None in the conversion function. It isn't that we don't know how to convert them, it is that we know we can't convert them.

Side note: You can remove the pass line of the except block's now.

@snowman2
Copy link
Member Author

snowman2 commented Apr 9, 2019

I like the idea of having a separate exclusions list in there for known skipped PROJ params. I think having this could also allow users to lookup which params are not supported by CF. It would probably also be good to have a separate one for known CF parameters not supported by PROJ. I think raising a warning for parameters not in either of those lists would be useful. But, open to opinions.

Another question in my mind is: Would it be better to only ignore parameters such as "nodefs" & "type" and warn for other ones such as "zone" or "k"?

You can remove the pass line of the except block's now.

That's neat. When was the feature introduced?

@snowman2
Copy link
Member Author

snowman2 commented Apr 9, 2019

Or, we could follow the current method in pyproj with an errcheck parameter to turn on/off warnings and then let the user turn on the warning to check and see what parameters are skipped with the default of being off. This would be a simpler route and would alleviate the need to manage more lists.

@djhoese
Copy link
Contributor

djhoese commented Apr 9, 2019

I'm interested in the idea of ignoring for some, warning for others. I could see it being a fine line though. I like the idea of an option for the user too.

For the pass thing, I only meant that you added skipped_params.append(proj_param) to the except: part of the code so you no longer need the pass which you added in your previous version of the code.

@snowman2
Copy link
Member Author

snowman2 commented Apr 9, 2019

I am leaning towards the errcheck option. Simple is always nice.

For the pass thing, ..

Ahhh, makes sense now. Thanks for the clarification :)

@djhoese
Copy link
Contributor

djhoese commented Apr 9, 2019

I was talking with @mraspaud about this and he pointed out the differences in oblique mercator handling. We had been using "alpha" instead of "gamma" for "azimuth_of_central_line". Also, we are using "lonc" for "longitude_of_projection_origin". Any idea if these can or should be used?

https://github.com/pytroll/satpy/blob/master/satpy/writers/cf_writer.py#L39-L55

@snowman2
Copy link
Member Author

snowman2 commented Apr 9, 2019

Thanks for bringing those points up. I am not entirely sure, so I will give my reasoning for my decisions and hopefully they will help choose which way to go on this. If you wouldn't mind sharing your reasoning, that would be helpful as well.

azimuth_of_central_line

Okay, for azimuth_of_central_line, it is defined in the Table F.1 (in the CF docs) as:

Specifies a horizontal angle measured in degrees clockwise from North. Used by certain projections (e.g., Oblique Mercator) to define the orientation of the map projection relative to a reference direction.

From the oblique mercator docs here: https://proj4.org/operations/projections/omerc.html

+gamma

 Azimuth of centerline clockwise from north of the rectified bearing of centre line. If +alpha is not given, then +gamma is used to determine +alpha.

+alpha

Azimuth of centerline clockwise from north at the center point of the line. If +gamma is not given then +alpha determines the value of +gamma.

longitude_of_projection_origin

For the longitude_of_projection_origin, I pulled lon_0 from Table 3 here: https://github.com/cf-convention/cf-conventions/wiki/Mapping-from-CF-Grid-Mapping-Attributes-to-CRS-WKT-Elements.

@snowman2
Copy link
Member Author

snowman2 commented Apr 9, 2019

Looks like y'all are using longitude_of_projection_origin=proj_dict.get('lon_0'), as well for another projection. Is it a special case for omerc to use lonc?

@snowman2
Copy link
Member Author

snowman2 commented Apr 9, 2019

Also looks like sweep=sweep_axis mapping is missing in the mappings in this PR.

@djhoese
Copy link
Contributor

djhoese commented Apr 9, 2019

According to https://proj4.org/operations/projections/omerc.html#cmdoption-arg-lonc lonc is the preferred reference longitude for omerc. And yes the sweep angle axis would be nice.

I'll let @mraspaud respond to the alpha versus gamma question. Not sure he has a ton of time though.

@snowman2
Copy link
Member Author

snowman2 commented Apr 9, 2019

From reading the docs you posted:

In the first case, the azimuth is given directly, using the option +alpha and defining the centre of projection using the options +lonc and +lat_0

It seems similar to the way the CF version is implementing it, so it sounds like the way you have it is the correct way.

@djhoese
Copy link
Contributor

djhoese commented Apr 9, 2019

@mraspaud said we wanted geotiff compatibility so he was looking at this page: http://geotiff.maptools.org/proj_list/hotine_oblique_mercator.html

@mraspaud
Copy link
Contributor

mraspaud commented Apr 9, 2019

@snowman2 sorry for coming in late on this one! So as @djhoese said, geotiff was the priority, and in the documentation they mention lonc and alpha. However, the CF document doesn't make it clear which one should be used. Can we see who pushed the omerc to the CF docs and ask about their thoughts on the matter ?

@snowman2
Copy link
Member Author

snowman2 commented Apr 9, 2019

@painter1, we are working on mapping the CF parameters to PROJ parameters and would appreciate it if you could clarify our questions about lonc and alpha with the oblique_mercator CF projection in this issue. Thanks!

@snowman2
Copy link
Member Author

Looking at the PROJ code here, it seems like alpha is the default and gamma is the backup. Also lonc is mentioned without lon_0 being referenced. So, I think it is likely that lonc and alpha are the ones to use.

@djhoese
Copy link
Contributor

djhoese commented Apr 10, 2019

Recent changes look good to me. Thanks.

@snowman2
Copy link
Member Author

Thanks for the review @djhoese, you had some good points/catches.

@snowman2
Copy link
Member Author

@jswhit are your concerns about warning the user resolved with the current implementation?

@jswhit
Copy link
Collaborator

jswhit commented Apr 11, 2019

@jswhit are your concerns about warning the user resolved with the current implementation?

Yes - that looks good to me.

@snowman2 snowman2 merged commit aa161c4 into pyproj4:master Apr 11, 2019
@snowman2
Copy link
Member Author

snowman2 commented May 2, 2019

Update: I am now very confident that lonc and alpha are the way to go with the coordinate_operation addition:

>>> from pyproj import CRS
>>> crs = CRS({
...         "proj": "omerc",
...         "lat_0": 10,
...         "lonc": 15,
...         "alpha": 0.35,
...         "gamma": 0.55,
...         "k": 1,
...         "x_0": 0,
...         "y_0": 0,
...         "ellps": "WGS84",
...         "units": "m",
...         "no_defs": None,
...         "type": "crs",
...     })
>>> crs.coordinate_operation
CONVERSION["unknown",
    METHOD["Hotine Oblique Mercator (variant B)",
        ID["EPSG",9815]],
    PARAMETER["Latitude of projection centre",10,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8811]],
    PARAMETER["Longitude of projection centre",15,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8812]],
    PARAMETER["Azimuth of initial line",0.35,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8813]],
    PARAMETER["Angle from Rectified to Skew Grid",0.55,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8814]],
    PARAMETER["Scale factor on initial line",1,
        SCALEUNIT["unity",1],
        ID["EPSG",8815]],
    PARAMETER["Easting at projection centre",0,
        LENGTHUNIT["metre",1],
        ID["EPSG",8816]],
    PARAMETER["Northing at projection centre",0,
        LENGTHUNIT["metre",1],
        ID["EPSG",8817]]]
>>> crs = CRS({
...             "proj": "omerc",
...             "lat_0": 10,
...             "lon_0": 15,
...             "alpha": 0.35,
...             "gamma": 0.35,
...             "k": 1,
...             "x_0": 0,
...             "y_0": 0,
...             "ellps": "WGS84",
...             "units": "m",
...             "no_defs": None,
...             "type": "crs",
...         })
>>> crs.coordinate_operation
CONVERSION["unknown",
    METHOD["Hotine Oblique Mercator (variant B)",
        ID["EPSG",9815]],
    PARAMETER["Latitude of projection centre",10,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8811]],
    PARAMETER["Longitude of projection centre",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8812]],
    PARAMETER["Azimuth of initial line",0.35,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8813]],
    PARAMETER["Angle from Rectified to Skew Grid",0.35,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8814]],
    PARAMETER["Scale factor on initial line",1,
        SCALEUNIT["unity",1],
        ID["EPSG",8815]],
    PARAMETER["Easting at projection centre",0,
        LENGTHUNIT["metre",1],
        ID["EPSG",8816]],
    PARAMETER["Northing at projection centre",0,
        LENGTHUNIT["metre",1],
        ID["EPSG",8817]]]

@snowman2 snowman2 deleted the cf branch September 2, 2019 14:47
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.

4 participants