[Proj] Pipeline +init files, inverse

[Thomas Knudsen]

... following up to myself here, since I seem to have gotten the inversions
wrong.

If I read your original init file correctly, what you want to do is this:

<case2>
proj=pipeline
step proj=vgridshift grids=OSGM15_GTX_OSGBtoGRS80.gtx
step proj=hgridshift grids=OSTN15_NTv2_OSGBtoETRS.gsb inv
step proj=tmerc lat_0=49 lon_0=-2 k=0.9996012717
                x_0=400000 y_0=-100000 ellps=airy

i.e. the hgridshift inverted, the other steps forward.

Saying
    echo -3 55 100 | cct -t0 +init=OSGBtrans:case2

Now results in:
      336128.9296    567726.4842       47.5270        0.0000



2018-02-23 11:00 GMT+01:00 Thomas Knudsen <knudsen.thomas at gmail.com>:

> Steve,
>
> Essentially (leaving out the double inversions), your "case 2" does this:
>
> <case2>
> proj=pipeline
> step proj=tmerc lat_0=49 lon_0=-2 k=0.9996012717
>                           x_0=400000 y_0=-100000 ellps=airy
> step proj=vgridshift grids=OSGM15_GTX_OSGBtoGRS80.gtx
> step proj=hgridshift grids=OSTN15_NTv2_OSGBtoETRS.gsb inv
>
> But your grid is given in geographic coordinates, not transverse mercator
> (cf. gdalinfo output [1] below), so starting your pipeline going from geo
> to transverse mercator takes you way outside the grid coverage.
>
> So start doing the gridshifts, then do the projection. When doing that, you
> get this:
>
> $ echo -3 55 100 | cct -t0 +init=OSGBtrans:case2
>   335948.6065    567743.2288       47.5270        0.0000
>
> /Thomas
>
>
>
>
> [1]  gdalinfo OSGM15_GTX_OSGBtoGRS80.gtx
>
> Driver: GTX/NOAA Vertical Datum .GTX
> Files: OSGM15_GTX_OSGBtoGRS80.gtx
> Size is 661, 1441
> Coordinate System is:
> GEOGCS["WGS 84",
>     DATUM["WGS_1984",
>         SPHEROID["WGS 84",6378137,298.257223563,
>             AUTHORITY["EPSG","7030"]],
>         AUTHORITY["EPSG","6326"]],
>     PRIMEM["Greenwich",0,
>         AUTHORITY["EPSG","8901"]],
>     UNIT["degree",0.0174532925199433,
>         AUTHORITY["EPSG","9122"]],
>     AUTHORITY["EPSG","4326"]]
> Origin = (-9.008333333333329,61.004166666666627)
> Pixel Size = (0.016666666666667,-0.008333333333333)
> Corner Coordinates:
> Upper Left  (  -9.0083333,  61.0041667) (  9d 0'30.00"W, 61d 0'15.00"N)
> Lower Left  (  -9.0083333,  48.9958333) (  9d 0'30.00"W, 48d59'45.00"N)
> Upper Right (   2.0083333,  61.0041667) (  2d 0'30.00"E, 61d 0'15.00"N)
> Lower Right (   2.0083333,  48.9958333) (  2d 0'30.00"E, 48d59'45.00"N)
> Center      (  -3.5000000,  55.0000000) (  3d30' 0.00"W, 55d 0' 0.00"N)
> Band 1 Block=661x1 Type=Float32, ColorInterp=Undefined
>   NoData Value=-88.8888000000000034
>
>
>
>
>
> 2018-02-23 8:46 GMT+01:00 stevejking <stevejking at talktalk.net>:
>
>> Kristian,
>> Thanks for the reply, I'll take a look at the link.
>>
>> For the moment, I've added the missing +step to case 3, and it now works
>> :)
>> I also removed the command line -t0's which I had in error put in the init
>> file.
>>
>> Unfortunately case 2 is still giving the error. I even tried removing the
>> breakline \'s.
>>
>> Will look more when time allows.
>>
>> ...
>>
>> Steve
>>
>>
>>
>> --
>> Sent from: http://osgeo-org.1560.x6.nabble.com/PROJ-4-f3840930.html
>> _______________________________________________
>> Proj mailing list
>> Proj at lists.maptools.org
>> http://lists.maptools.org/mailman/listinfo/proj
>>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.maptools.org/pipermail/proj/attachments/20180223/e732de6a/attachment.htm