2

We are using pyproj to transform between coordinate reference systems. What is a good way to check if the transformation is simply scale/translate/units? Basically, I want to know if a grid in the first crs will still be in a grid in the target crs.

Edit: I want to add some clarity about what I am trying to detect. Suppose I have a 3x4 grid of 12 points represented by 3 x values and 4 y values. When I transform these points, sometimes the grid is distorted such that 12 x,y pairs are required. But sometimes, the transformation yields a 3x4 grid that can still be represented using 3 x values and 4 y values. It is these latter cases that we would like to detect.

An example is EPSG:4236 to EPSG:3395. Other transforms distort the grid, such as EPSG:4236 to EPSG:2193.

I've listed out a lot of detail below, but I think this is going to be a pretty easy question. Ideally, there is something simple in the pyproj.Transformer object that I can inspect, or an easy way to compare projection definitions in some format.

  1. Empirical Test

    Obviously, I can transform a few points and then check if they are still gridded. I'd like to find a better way.

  2. Compare DATUM

    Using pyproj.CRS('EPSG:4236').datum, the datums are

    • EPSG:4236: DATUM["World Geodetic System 1984", ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]], ID["EPSG",6326]]
    • EPSG:3395: DATUM["World Geodetic System 1984", ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]], ID["EPSG",6326]]
    • EPSG:2193: DATUM["New Zealand Geodetic Datum 2000", ELLIPSOID["GRS 1980",6378137,298.257222101, LENGTHUNIT["metre",1]], ID["EPSG",6167]]

    So the DATUM for 4326 and 3395 are the exact same. Is it the case that two crs with the same exact datum string is always a scale/translate? Is it enough to compare just the DATUM portion (World Geodetic System 1984 vs New Zealand Geodetic Datum 2000) or just the ID portion (EPSG 6326 vs EPSG 6167)? Or, is this not going to work at all (can two CRS share the same DATUM but still require distortion when transforming)?

    Furthermore, even if that does work, does every crs have a "datum"? It would be nice to detect this even with custom crs.

  3. Compare PROJ4

    The proj4 strings:

    • EPSG:4326: +proj=longlat +datum=WGS84 +no_defs +type=crs
    • EPSG:3395: +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs
    • EPSG:2193: +proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +units=m +no_defs +type=crs

    Anything I can do with this? Seems like there is not quite sufficient info, here, but maybe some values have "defaults"? Which parts of these strings indicate scaling/transformation as opposed to other computations?

  4. Compare WKT

    The WKT strings or JSON representations have the same basic parameters, but they are more detailed and may be simpler to use. I won't post them here entirely. But, for example, EPSG:3395 has EPSG:4326 as the base_crs, but of course EPSG:2193 does not. Would comparing base_crs be sufficient?

  5. Inspect the pyproj Transformer object

    I'm using a pyproj Transformer, so maybe the best thing to do is simply inspect it.

    • Using the transformer definition directly might be better than trying to compare proj4 strings; which part(s) should I use?
    • Alternatively, the operations is empty for the 4236 to 3395 transformation, but not empty for 4326 to 2193. Is it as simple as that?
    >>> t = pyproj.Transformer.from_crs('EPSG:4236', 'EPSG:3395')
    >>> t.definition
    'proj=pipeline step proj=axisswap order=2,1 step proj=unitconvert xy_in=deg xy_out=rad step proj=merc lon_0=0 k=1 x_0=0 y_0=0 ellps=WGS84'
    >>> t.operations
    ()
    
    >>> t = pyproj.Transformer.from_crs('EPSG:4236', 'EPSG:2193')
    >>> t.definition
    'proj=pipeline step proj=axisswap order=2,1 step proj=unitconvert xy_in=deg xy_out=rad step proj=tmerc lat_0=0 lon_0=173 k=0.9996 x_0=1600000 y_0=10000000 ellps=GRS80 step proj=axisswap order=2,1'
    >>> t.operations
    (<Coordinate Operation: Transformation>
     Inverse of NZGD2000 to WGS 84 (1)
     Area of Use:
     - name: New Zealand
     - bounds: (160.6, -55.95, -171.2, -25.88), <Coordinate Operation: Conversion>
     New Zealand Transverse Mercator 2000
     Area of Use:
     - name: New Zealand - onshore
     - bounds: (166.37, -47.33, 178.63, -34.1))
    

1 Answer 1

2

As a quick answer, most transforms are more complex than a scale/translate/unit, so you are looking for exceptions. The distortion come from the change in datum (that is the transformation, from lat/long to lat/long) and from the type of projection (from lat/long to XY).

Having the same (or rotated around the vertical) datum is a prerequisite but, if you have different projections on top of the datum, this will "distort" your regular grid.

Edit: as a remark, your square in EPSG4326 will become rectangles in epsg3395 (elongation depending of the latitude). If you don't mind those different rescaling in X and Y directions, then all the cylindrical projections with a vertical axis will yield a grid of rectangles (squares in the special case of the Plate Carrée). Unfortunately the name "cylindrical" is not always in the description of the CRS, so you would need to build a look up table first).

9
  • Yep, I'm looking for the exceptions. And yep, rectangles are okay, as long as it is still a grid. In fact, a rotated grid is okay, too.
    – jmilloy
    Feb 4, 2020 at 16:48
  • So, do you know how to tell if the projection on top of the datum "distorts" the grid or just scales/translates it?
    – jmilloy
    Feb 4, 2020 at 16:52
  • 1
    The 'grid' in 3395 scales differently at different latitudes.If you're working in a small area, then single scale values for the x and y coordinates might work. That would be generally true for many PCS--has to be a small area. I have a feeling that the code would spend as much time calculating the scales/rotations/translations / identifying the exceptions than just performing the reprojections.
    – mkennedy
    Feb 5, 2020 at 0:39
  • @mkennedy I can't produce what you are claiming using pyproj. It seems to be scaling the same at every latitude.
    – jmilloy
    Feb 5, 2020 at 19:08
  • @radouxju The lookup table doesn't seem the right solution for me.
    – jmilloy
    Feb 5, 2020 at 19:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.