Use input time in epsg transform?

[lucasw]

With pyproj Transformer I had this:

transformer = Transformer.from_crs(6319, 7912, always_xy=True, allow_ballpark=False)

and then providing the year (with a fractional component in tt_year:

 x, y, z, tt = self.transformer.transform(xx=lon, yy=lat, zz=alt, tt=tt_year)

In rust

let from = "EPSG:6319";
let to = "EPSG:7912";
let proj_convert = proj::Proj::new_known_crs(from, to, None).unwrap();
...
let lon_lat = proj_convert.convert((lon, lat)).unwrap();

I think it's using 'right now' internally as the year time, if I make the python use right now the output seems very close. But can I craft a proj string or do something else to be able to input the time in rust? (matching the allow_ballpark would be desirable as well)

[urschrei]

It looks like the ballpark setting is one of the new supported OPTIONS: https://proj.org/en/stable/development/reference/functions.html#c.proj_create_crs_to_crs_from_pj

I'm less clear on the timestamp – is it part of the CoordinateMetadata discussed in the link above? The 'velocity grids' note makes me doubt that, though. In any case: we should probably support the OPTIONS setting fully, and this looks fairly straightforward, but I'm not sure about the timestamp issue – that's almost certainly a larger piece of work requiring some thought and discussion around design and implementation (would this require a datetime library as a new dependency? What would we use? etc), and we need some more information first: is this related to CoordinateMetadata?

[urschrei]

@lucasw #226 and #227 have some progress on this: #226 allows you to add epoch metadata (e.g. 2010.0) to an existing transform using Proj::coordinate_metadata_create, and #227 allows you to use the resulting Proj instance directly as input to a new transformation pipeline using Proj::create_crs_to_crs_from_pj (see the PR source for examples etc). Without #227 I think you can still use the epoch-enabled objects by e.g. by using their string representations (via Proj::def()) as input to to Proj::new_known_crs()

Pyproj uses a slightly different API than the libproj API, so I think I was a little confused, but there seems to be only one canonical function in libproj to explicitly add an epoch (coordinate_metadata_create), which we're now providing access to. Please take a look and let me know if this works.

[lucasw]

Using #226 test_metadata_creation() if I set this

        let wgs84 = "EPSG:6319";
        let to = "EPSG:7912";
        let proj = Proj::new_known_crs(wgs84, to, None).unwrap();

        let epoch = 2021.3;
        let np = proj.coordinate_metadata_create(epoch).unwrap();
        assert_eq!(np.coordinate_metadata_get_epoch(), epoch);

I get

called `Result::unwrap()` on an `Err` value: ProjError("Unknown error (code 4096)")

for that coordinate_metadata_create() line

[urschrei]

Hmm I need to check whether libproj also does this, but the problem with the above is that our API is currently too general: it allows calling coordinate_metadata_create on Proj instances which are definitions (what you get if you call Proj::new()) AND transformations (what you get if you call Proj::new_crs(). But coordinate_metadata_create only makes sense for the former: I don't think a transformation pipeline can have an epoch. Only a transformation object. This passes on #227 (would be good to test it with some coordinates and a valid epoch and see whether results differ):

    #[test]
    fn metadata_bad() {
        let wgs84 = Proj::new("EPSG:6319").unwrap();
        let to = Proj::new("EPSG:7912").unwrap();
        let epoch = 2021.3;
        let with_epoch = to.coordinate_metadata_create(epoch).unwrap();

        let proj_without_epoch = wgs84.create_crs_to_crs_from_pj(&to, None, None).unwrap();
        let proj_epoch = wgs84
            .create_crs_to_crs_from_pj(&with_epoch, None, None)
            .unwrap();

        assert_ne!(proj_without_epoch.def().unwrap(), proj_epoch.def().unwrap());
    }

Anyway, I need to see whether libproj also allows the ability to "hold it wrong", and to figure out whether we can add this in a way that prevents that, without a breaking API change.

[lucasw]

Ok using the #227 code and the metadata_bad() example I'm able to get within 1e-9 of the numbers I was expecting (vs 1e-6 degrees with no epoch specified), which is about what python was doing as well, and it seems good enough- though I'm curious why python and rust/c++ don't match exactly (different proj versions?).

[urschrei]

It may be different libproj versions and / or the presence of correction grids.

To ensure that the transformation pipeline is using grids when they're available for a given transform you can use the network feature and then activate it by calling enable_network(true) on a ProjBuilder instance (see https://docs.rs/proj/latest/proj/#network-cache-and-search-path-functionality).

Something like (pull #227 again first to get the ProjBuilder update):

let wgs84 = Proj::new("EPSG:6319").unwrap();
let to = Proj::new("EPSG:7912").unwrap();
let epoch = 2021.3;
let with_epoch = to.coordinate_metadata_create(epoch).unwrap();


let mut builder = ProjBuilder::new();
builder.enable_network(true).unwrap();
// Disable caching to ensure we're accessing the network.
builder.grid_cache_enable(false);

let transformer = builder
    .proj_create_crs_to_crs_from_pj(&from, &with_epoch, None, None)
    .unwrap();
let result = transformer
    .convert(MyPoint::new(4760096.421921, 3744293.729449))
    .unwrap();

I don't know whether there's a grid available for that transformation pair, but worth trying.

[urschrei]

This functionality is now available in v0.30.0, released today, so I'm going to close this.