2

I am trying to transform coordinates from a local grid to a projected coordinate system (UTM) to ultimately reproject raster or vector data.

I am using :

  • Python 3.9
  • GDAL 3.6.0
  • Pyproj 3.4

I generated a test point dataset with local and real world coordinates in NAD83 UTM Z13 (EPSG:26913) for each point.

Three points from the dataset have been used to compute affine parameters to define a WKT2 transform based on the posts below:

how-to-use-affine-transform-parameters-in-a-proj-or-wkt2-string

how-do-i-find-these-values-for-an-affine-transformation-state-plane-to-local-gr

Base coordinates in NAD83 UTM Z13 (EPSG:26913)

    base = [[468000,6463900,0],[483000,6481880.76,0],[496856.41,6481880.76,0]]

Local grid coordinates

    local = [[0,0,0],[30000,0,0],[30000,16000,0]]

And the WKT2_2018 format affine transformation obtained from these points is:

    DERIVEDPROJCRS["survey_grid",BASEPROJCRS["NAD83 / UTM zone 13N",BASEGEOGCRS["NAD83",
    DATUM["North American Datum 1983",ELLIPSOID["GRS 1980",6378137,298.257222101,LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4269]],
    CONVERSION["UTM zone 13N",METHOD["Transverse Mercator",ID["EPSG",9807]],
    PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],
    PARAMETER["Longitude of natural origin",-105,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],
    PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1],ID["EPSG",8805]],
    PARAMETER["False easting",500000,LENGTHUNIT["metre",1],ID["EPSG",8806]],
    PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]]],
    DERIVINGCONVERSION["Affine",METHOD["Affine parametric transformation",ID["EPSG",9624]],
    PARAMETER["A0",-6259712.633486367,LENGTHUNIT["metre",1],ID["EPSG",8623]],
    PARAMETER["A1",0.5421061154648178,SCALEUNIT["coefficient",1],ID["EPSG",8624]],
    PARAMETER["A2",0.929306188894481,SCALEUNIT["coefficient",1],ID["EPSG",8625]],
    PARAMETER["B0",3067450.732242991,LENGTHUNIT["metre",1],ID["EPSG",8639]],
    PARAMETER["B1",0.929306188894481,SCALEUNIT["coefficient",1],ID["EPSG",8640]],
    PARAMETER["B2",-0.5421061154648178,SCALEUNIT["coefficient",1],ID["EPSG",8641]]],
    CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",north,
    ORDER[2],LENGTHUNIT["metre",1]]]

I am now trying to use this WKT transformation pipeline with Python to reproject raster or vector data using the GDAL Python API but I can't seem to understand how to implement it.

The GDAL C API function OCTNewCoordinateTransformationEx creates a new transformation object that can use a transform pipeline but I don't see its equivalent in the Python API.

If I was to consider points only, I could use Pyproj v3.4 as it can define a transform from a pipeline using Transformer.from_pipeline.

Any recommendation on how to implement this with the GDAL Python API?

user1639843
  • 81
  • 5
  • I think you need 3 match-up points to get all 6 affine parameters. – Pointdump Dec 14 '22 at 23:56
  • The question is more on how to use the affine parameters and not about the robustness of the transformation. I used the minimum number of points required by simil.py but 3 points would indeed provide better results ! I upated the post adding a point and updating the affine parameters. – user1639843 Dec 20 '22 at 01:10
  • 1
    With gdalwarp and ogr2ogr the affine transformation can be used in a +proj=affine pipeline with the -ct parameter. If you want to transform from the local to the projected (utm 13N) srs, the transformation must be the inverse of the one used in the derived from projected wkt. And the target srs must be defined with a -t_srs parameter. I don't know if you can use the derived wkt as source srs and save the pipeline, i think for compatibility reasons gdal still needs wkt1 srs. I haven't made a Python script with that process yet, in order to answer accordingly. – Gabriel De Luca Dec 20 '22 at 11:07

0 Answers0