TransWikia.com

Reproject WGS84 raster file to topocentric LTP-ENU using GDAL with a custom proj string

Geographic Information Systems Asked on February 15, 2021

I have some raster files which are given in WGS84. As an example, here are some gdalinfo for one of those:

Driver: GTiff/GeoTIFF
Files: dem_r.tif
Size is 2000, 911
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 = (-43.795421989932834,-22.745462070999945)
Pixel Size = (0.000371130000000,-0.000371130000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -43.7954220, -22.7454621) ( 43d47'43.52"W, 22d44'43.66"S)
Lower Left  ( -43.7954220, -23.0835615) ( 43d47'43.52"W, 23d 5' 0.82"S)
Upper Right ( -43.0531620, -22.7454621) ( 43d 3'11.38"W, 22d44'43.66"S)
Lower Right ( -43.0531620, -23.0835615) ( 43d 3'11.38"W, 23d 5' 0.82"S)
Center      ( -43.4242920, -22.9145118) ( 43d25'27.45"W, 22d54'52.24"S)
Band 1 Block=2000x1 Type=Float32, ColorInterp=Gray
  NoData Value=-3.4028234663852886e+38

I am now trying to reproject it to a local tangent plane oriented East, North, Up (LTP-ENU), also known as a topocentric coordinate reference system. To this end I naturally went to gdalwarp where I discovered this relatively new -ct flag which seems promising:
The GDAL -ct flag

Therefore, I constructed this proj string as follow:

projstring="+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m
  +step +proj=cart +ellps=WGS84
  +step +proj=topocentric +lat_0=-22.7454621 +lon_0=-43.7954220 +h_0=0 +ellps=WGS84
"

When I ran gdalwarp, everything went fine (I hope):

$ gdalwarp -ct "${projstring}" dem_r.tif dem_r_warped.tif
Creating output file that is 975P x 1975L.
Processing dem_r.tif [1/1] : 0Using internal nodata values (e.g. -3.40282e+38) for image dem_r.tif.
Copying nodata values from source dem_r.tif to destination dem_r_warped.tif.
...10...20...30...40...50...60...70...80...90...100 - done.

But then gdalinfo gives me:

$ gdalinfo dem_r_warped.tif
Driver: GTiff/GeoTIFF
Files: dem_r_warped.tif
Size is 975, 1975
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 = (-37632.325666819946491,76238.183460185828153)
Pixel Size = (38.608857984279830,-38.608857984279830)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  -37632.326,   76238.183) (Invalid angle,Invalid angle)
Lower Left  (  -37632.326,     -14.311) (Invalid angle, 14d18'39.81"S)
Upper Right (      11.311,   76238.183) ( 11d18'39.12"E,Invalid angle)
Lower Right (  11.3108679, -14.3110588) ( 11d18'39.12"E, 14d18'39.81"S)
Center      (  -18810.507,   38111.936) (Invalid angle,Invalid angle)
Band 1 Block=975x2 Type=Float32, ColorInterp=Gray
  NoData Value=-3.4028234663852886e+38

Two details struck me:

  1. The size seems reversed Size is 975, 1975. It is also confirmed when I open the raster in QGIS; it looks like it has first been mirrored horizontally along a vertical axis passing by the middle of the raster, and then rotated 90° clockwise.
  2. Some corner coordinates gives (Invalid angle,Invalid angle).

Maybe all these has to do with correctly formatting the proj string, but then how to do it correctly and how does it deal with negative values for latitude and longitude (because I have rasters all over the world, in the northern and the southern hemisphere)? I’m lost for a while in the proj doc without any success and a trial and error approach doesn’t help me much in understanding what’s happening under the hood.

What I was expecting is something like:

$ gdalinfo dem_r_warped.tif
Driver: GTiff/GeoTIFF
Files: dem_r_warped.tif
Size is 975, 1975
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 = (0,0,0)
Pixel Size = (1,1)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -38119.890905  18672.743845 -141.335023)
Lower Left  ( -38025.461332 -18769.163170 -141.054109)
Upper Right (  38119.890905  18672.743845 -141.335023)
Lower Right (  38025.461332 -18769.163170 -141.054109)
Center      (             0             0           0)
Band 1 Block=975x2 Type=Float32, ColorInterp=Gray
  NoData Value=-3.4028234663852886e+38

Note that here, all corners were manually computed relative to the center with:
echo <easting_corder>E <northing_corner>N 0 | CartConvert -l <easting_center>E <northing_center>N 0

So, origin is (0,0,0) and pixel size is (1,1) (this latter value must be the result of a resampling step).

Here are some environment information:

# gdalinfo --version
GDAL 3.3.0dev-28306791ede238f0076ae83c919d245a92907d6f, released 2021/01/15
# proj --version
Rel. 8.0.0, March 1st, 2021

Ubuntu 20.04.1 LTS:

# uname -a  
Linux a5de08b82589 5.4.0-60-generic #67~18.04.1-Ubuntu SMP Tue Jan 5 22:01:05 UTC 2021 x86_64 x86_64 x86_64 GNU/Linux

This thread also helped me a lot up to that point.

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP