[PROJ] Affine transformation in proj < 6?

[Tammo Jan Dijkema]

Thanks for the replies! I realised later on that the Earth is locally quite flat, so I described in the geotiff only the latitude and longitude of the lower left corner. The metadata doesn't need to know everything about my coordinate system.

So I did the conversion from my coordinate system to lat-long before writing anything to tiff. There is an affine transform stored in geotiff.

I found the python library rasterio, which I use as follows. This does the trick for me:

    urc = pqr_to_latlongheight(urc_pqr) # upper right corner
    llc = pqr_to_latlongheight(llc_pqr) # lower left corner

    transform = Affine.translation(llc[0] - long_res / 2,
                                   llc[1] - lat_res / 2) * \
                  Affine.scale(long_res, lat_res)

    with rasterio.open(filename, "w", driver="GTiff",
                       height=height, width=width,
                       count=1,
                       dtype=image.dtype,
                       crs='+proj=latlong',
                       transform=transform) as gtif:
        gtif.write(image, 1)

Note that there is a proj string in the GeoTiff, but a quite simple one. I wouldn't be surprised if without the proj string, it would work as well (not tested).

The file produced this way loads perfectly into QGIS.

Tammo Jan


> Op 12 dec. 2019, om 17:32 heeft Even Rouault <even.rouault at spatialys.com> het volgende geschreven:
> 
> On jeudi 12 décembre 2019 13:31:37 CET Tammo Jan Dijkema wrote:
>> Hi all,
>> 
>> I'm new to proj, so sorry if this is a trivial question (I've read through
>> docs for an hour without finding an answer, so it can't be too trivial).
>> 
>> I have a custom coordinate system, which is defined through an affine
>> transformation from ETRS89. All coordinates I care about are in xyz
>> Cartesian coordinates. In particular, to get pqr coordinates, I do
>> 
>> coords_pqr = (coords_etrs - pqrorigin_etrs) * (3x3 matrix)
>> 
>> Now I'm generating raster images in PQR, which I would like to store in
>> GeoTiff, so if I understand correctly I need to write a proj string into
>> the GeoTiff to describe the projection.
> 
> Except that GeoTIFF doesn't store PROJ strings...
> Provided that you manage a PROJ string to expression what you want, you could 
> (ab)use a GDAL WKT1 string to embed a PROJ string in it, and embed that WKT 
> string into GeoTIFF, but the compatibility of this will be really limited.
> 
> Example of this:
> $ gdal_translate in.tif foo.tif -a_srs "+proj=geos +sweep=x +datum=WGS84 
> +h=123456789"
> 
> $ listgeo foo.tif
> 
>      PCSCitationGeoKey (Ascii,630): "ESRI PE String = 
> PROJCS["unknown",GEOGCS["unknown",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"]]],PROJECTION["Geostationary_Satellite"],PARAMETER["central_meridian",
> 0],PARAMETER["satellite_height",123456789],PARAMETER["false_easting",
> 0],PARAMETER["false_northing",0],UNIT["metre",
> 1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],EXTENSION["PROJ4","+proj=geos 
> +sweep=x +lon_0=0 +h=123456789 +x_0=0 +y_0=0 +datum=WGS84 +units=m 
> +no_defs"]]"
> 
> What you could do instead of abusing the CRS is to provide a geotransform 
> matrix with the affine transformation (assuming it is only in the X,Y 
> projected coordinate space)