openmethane / openmethane-prior

Method to calculate a gridded, prior emissions estimate for methane across Australia
Apache License 2.0
4 stars 0 forks source link

Domain projection may be missing x/y offset/origin #24

Open aethr opened 1 month ago

aethr commented 1 month ago

Describe the bug

When reading the domain file in omInputs.py, a pyproj Proj object is instantiated using the parameters from the domain netCDF file. However, the grid offset was only recently added to the domain file, and is not used when instantiating the projection.

        domainProj = pyproj.Proj(
            proj='lcc',
            lat_1=domainXr.TRUELAT1,
            lat_2=domainXr.TRUELAT2,
            lat_0=domainXr.MOAD_CEN_LAT,
            lon_0=domainXr.STAND_LON,
            a=6370000,
            b=6370000,
        )

This is fine when using the projection directly, but when using it as a CRS this may cause projected coordinates to be in totally the wrong place.

The pyproj docs describe x_0 and y_0 parameters to a lambert projection: https://proj.org/en/9.4/operations/projections/lcc.html#lcc

I'm not 100% clear on the terms, but I believe this may map to the XORIG and YORIG attributes in the domain file.

Failing Test

I'm not 100% sure about the math, but in theory if we add the x/y coordinate origin to the projection to make a CRS, and then use it to transform the same numbers, that should take us back to the grid centre-point?

def test_crs_coordinates():
    toLatLong = Transformer.from_crs(domainProj.crs, 4326)

    centreLatLong = list(toLatLong.itransform(
        points=[(domainXr.attrs['XORIG'], domainXr.attrs['YORIG'])],
        switch=True,
    ))[0]

    assert centreLatLong == (domainXr.attrs['XCENT'], domainXr.attrs['YCENT'])

-->

Expected behavior

Using domainProj.crs for transformations should result in lat / long values in the right position on the earth w.r.t. our grid points.

Screenshots

System

Additional context

May be the cause of #17 .

aethr commented 6 days ago

I can confirm that naively adding x_0 and y_0 when constructing the pyproj projection in https://github.com/openmethane/openmethane-prior/blob/main/src/openmethane_prior/omInputs.py#L89-L91 did not resolve the issue. This majorly distorted the attribution of emissions to grid cells in several layers.

prayner commented 6 days ago

The problem was also mentioned in issue #17. I've established it doesn't happen for all sectors but not yet found the common factor among the sectors where it does.