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 unit_conversion_factor for units in coordinate axis to represent the number of SI standard units per unit (i.e. meters). #248

Closed
snowman2 opened this issue Feb 14, 2020 · 12 comments
Labels
enhancement Proposals to add new capabilities, improve existing ones in the conventions, improve style or format

Comments

@snowman2
Copy link
Contributor

snowman2 commented Feb 14, 2020

Title

Add unit_conversion_factor for units in coordinate axis to represent the number of SI standard units per unit (i.e. meters).

Moderator

Any takers?

Moderator Status Review [last updated: YY/MM/DD]

???

Requirement Summary

http://docs.opengeospatial.org/is/12-063r5/12-063r5.html#35
Currently, in the WKT form, you can specify the conversion factor to meters along with the name of the units for a cartesian coordinate system:

CS[Cartesian,2],
    AXIS["easting (X)",east,
        ORDER[1],
        LENGTHUNIT["US survey foot",0.304800609601219,
            ID["EPSG",9003]]],
    AXIS["northing (Y)",north,
        ORDER[2],
        LENGTHUNIT["US survey foot",0.304800609601219,
            ID["EPSG",9003]]]

This is also the same for degrees to convert to radians in an ellipsoidal projection system:

CS[ellipsoidal,2],
    AXIS["geodetic latitude (Lat)",north,
        ORDER[1],
        ANGLEUNIT["degree",0.0174532925199433,
            ID["EPSG",9122]]],
    AXIS["geodetic longitude (Lon)",east,
        ORDER[2],
        ANGLEUNIT["degree",0.0174532925199433,
            ID["EPSG",9122]]]

Technical Proposal Summary

For example:

  double y(y);
    y:units = "US survey foot";
    y:unit_conversion_factor: 0.304800609601219;
    y:long_name = "y coordinate of projection";
    y:standard_name = "projection_y_coordinate";
  double x(x);
    x:units = "US survey foot";
    x:unit_conversion_factor: 0.304800609601219;
    x:long_name = "x coordinate of projection";
    x:standard_name = "projection_x_coordinate";

Benefits

This removes the need of having a lookup table to convert any unit the provider gives to meters (or whatever the standard unit is for the axis). It makes the provider of the dataset responsible for giving that information to the user.

Status Quo

I couldn't find something like this in the spec.

Detailed Proposal

See above ^^

@snowman2 snowman2 added the enhancement Proposals to add new capabilities, improve existing ones in the conventions, improve style or format label Feb 14, 2020
@zklaus
Copy link

zklaus commented Feb 14, 2020

Units have to come from udunits (see CF-1.8, 3.1), which will also provide the conversion factors. In this sense, this proposal seems unnecessary.

@martinjuckes
Copy link
Contributor

As Klaus points out, the units need to come from the Udunits list, so the correct units attribute would be US_survey_foot.

When I use the Udunits2 linux command line tool, it gives a conversion factor of 0.304801, which is significantly less precise. The XML database that comes with the python package gives the more precise value of 1200/3937.

Given that it takes some digging to get the conversion factor out of Udunits, I can see a potential case for providing space to add it explicitly. You would need, however, to indicate what the target units are, e.g. meter_conversion_factor.

@marqh
Copy link
Member

marqh commented Feb 14, 2020

Hello @snowman2

thank you for highlighting this detail clearly and comprehensively.

Given CF's dependent relationship with UDUnits, which is not specified by CF, I am tempted to concur with the position put forward by @zklaus (#248 (comment))

I think that not including unit conversion factors within CF metadata seems a proportionate and sensible position.
This approach would imply (to me) a separation of the WKT-CRS units specification and implementation from the UDUnits implementation that CF depends upon.

I am open to influence, but I don't see a compelling requirement here.

mark

@davidhassell
Copy link
Contributor

Hello, I also think that this is not a desirable change, since it creates a redundancy in information, which increases the likelihood of inconsistencies.

David

@snowman2
Copy link
Contributor Author

Another application where this would be useful is the geostationary projection.

How the conversion is handled for geostationary:

Related to:

Adding the conversion_factor to the axis information would remove the need for the software application or user to do have any special cases for the projection as it will just work if this spec was added.

For example, the solution could instead be resolved with:

import xarray
import rioxarray
from pyproj import CRS

xds = xarray.open_dataset("OR_ABI-L2-LSTF-M6_G17_s20200341900321_e20200341909388_c20200341910038.nc")
sat_height = xds.goes_imager_projection.attrs["perspective_point_height"]
crs_meter = CRS.from_cf(xds.goes_imager_projection.attrs)
goes_json = crs_meter.to_json_dict()

unit_deg = {"type": "LinearUnit", "name": "Deg", "conversion_factor": sat_height}
goes_json["coordinate_system"] = {
    "subtype": "Cartesian",
    "axis": [
      {
        "name": "Easting",
        "abbreviation": "E",
        "direction": "east",
        "unit": unit_deg
      },
      {
        "name": "Northing",
        "abbreviation": "N",
        "direction": "north",
        "unit": unit_deg
      }
    ]
}

cc = CRS(goes_json)
xds.rio.write_crs(cc, inplace=True)
xds = xds[["LST", "DQF"]]
xds3857 = xds.rio.reproject("EPSG:3857")
xds3857.rio.to_raster("geos17_epsg3857.tif")

I verified that this works:
image

This is not something that can be solved with a lookup table since the value varies depending on the satellite height.

@JonathanGregory
Copy link
Contributor

Like others, I would argue against this change, on grounds of wanting to avoid redundancy and hence inconsistency. A main purpose of providing the CRS information (in grid_mapping attributes or as WKT) is to enable the software to convert from the grid coordinates to other coordinate systems. It is expected that software which processes the data will be able to do that.

@JimBiardCics
Copy link
Contributor

I agree that we shouldn't make this change. It looks like a way to allow people to get extremely sloppy with units.

@snowman2
Copy link
Contributor Author

snowman2 commented Mar 2, 2020

I have a list of units I would like to ensure are supported at a minimum here to ensure a safe conversion from WKT to the CF convention (note: only the ones with the type length). See: unit_of_measure.sql.

And from what I can tell, only a subset are supported by the CF convention: UDUNITS Database.

Here are some examples of missing units:

  • British foot (Sears 1922) - 3.04799471538676203241e-01
  • Indian foot (1975) - 0.3047995
  • Gold Coast foot - 3.04799710181508809458e-01

I am sure there are others not supported by CF or in the EPSG databsse currently (WKT can support any unit as long as it has the unit conversion factor). When a unit that is not supported by the CF convention is found, what is the recommended procedure? Is this something that the developers of the UDUNITS program are interested in keeping in-sync?

@snowman2 snowman2 changed the title Add unit_conversion_factor for units in coordinate axis to convert to meters Add unit_conversion_factor for units in coordinate axis to represent the number of SI standard units per unit (i.e. meters). Mar 2, 2020
@semmerson
Copy link

@snowman2,

I'm the UDUNITS developer.

I've never heard of these units. Do you have a reference for them (other than the GitHub repository)?

@snowman2
Copy link
Contributor Author

snowman2 commented Mar 2, 2020

@semmerson, you can go to http://www.epsg-registry.org/ and filter by the unit of measure type. For example:
image

@JimBiardCics
Copy link
Contributor

@snowman2 It seems pretty clear that this is a UDUNITS issue, not a CF issue.

@snowman2
Copy link
Contributor Author

snowman2 commented Mar 2, 2020

I think that is all I have to offer for clarifying information on this issue. I will close it now since it seems that 👎 is the consensus here on the proposal.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement Proposals to add new capabilities, improve existing ones in the conventions, improve style or format
Projects
None yet
Development

No branches or pull requests

8 participants