[gdal-dev] Move pixel ULx ULy - what is dx_dy + dy_dx

[Andrew C Aitchison]

On Fri, 16 Dec 2022, Clive Swan wrote:

> Greetings,
>
> I want to shift the UL x, UL y of an image.
>
> I am trying to find out what the *dx_dy + dy_dx* is.

They represent (part of) the rotation;
technically I think that they are they are the "skew".
For a *translation* they are both equal to zero.

> Is dx_dy, dy_dx = (0,0)
> The Top Left/Right value??

For your transform, yes they are zero,
but no, not the Top Left/Right.

> The code is:
>
> import gdal
> ### open dataset with update permission
> ds = gdal.Open('/data/coastal-2020.tif', gdal.GA_Update)
> ### get the geotransform as a tuple of 6
> gt = ds.GetGeoTransform()
> ### unpack geotransform into variables
> x_tl = (-180.0000000)
> x_res = 0.000800000000000
> *dx_dy* =
> y_tl =  (90.0000000)
> *dy_dx* =
> y_res = -0.000800000000000
>
> x_tl, x_res, *dx_dy*, y_tl, *dy_dx*, y_res = gt
>
> # compute shift of 1 pixel RIGHT in X direction
> shift_x = 1 * x_res
> # compute shift of 2 pixels UP in Y direction
> # y_res likely negative, because Y decreases with increasing Y index
> shift_y = -2 * y_res
>
> # make new geotransform
> gt_update = (x_tl + shift_x, x_res, dx_dy, y_tl + shift_y, dy_dx, y_res)

x and y are in *degrees*.
You *can* shift your map left/right = east/west, but odd things may happen
at 180 degree longitude, since a common convention is that this
is the left-right border.
However, you cannot shift your map up/down = north/south as this will move
(in your case) the top of the map north of the north pole.
This is particularly bad with any Mercator projection,
since the north pole is already at infinity on the map.

>> From gdalinfo:
> Size is 450000, 225000
> Coordinate System is:
> GEOGCRS["WGS 84",
>    DATUM["World Geodetic System 1984",
>        ELLIPSOID["WGS 84",6378137,298.257223563,
>            LENGTHUNIT["metre",1]]],
>    PRIMEM["Greenwich",0,
>        ANGLEUNIT["degree",0.0174532925199433]],
>    CS[ellipsoidal,2],
>        AXIS["geodetic latitude (Lat)",north,
>            ORDER[1],
>            ANGLEUNIT["degree",0.0174532925199433]],
>        AXIS["geodetic longitude (Lon)",east,
>            ORDER[2],
>            ANGLEUNIT["degree",0.0174532925199433]],
>    ID["EPSG",4326]]
> Data axis to CRS axis mapping: 2,1
> Origin = (-180.000000000000000,90.000000000000000)
> Pixel Size = (0.000800000000000,-0.000800000000000)
>
> Image Structure Metadata:
>  COMPRESSION=LZW
>  INTERLEAVE=BAND
>  PREDICTOR=3
> Corner Coordinates:
> Upper Left  (-180.0000000,  90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"N)
> Lower Left  (-180.0000000, -90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"S)
> Upper Right ( 180.0000000,  90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"N)
> Lower Right ( 180.0000000, -90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"S)
> Center      (   0.0000000,   0.0000000) (  0d 0' 0.01"E,  0d 0' 0.01"N)

-- 
Andrew C. Aitchison                      Kendal, UK
                    andrew at aitchison.me.uk