Wrong contours when generating from geotiff with other projections

[rmanhaeve]

Hi

I tried generating contours from a geotiff file that had a Lambert2008 projection. They are wrong however, and do not match up with the ones (correctly) generated from the same geotiff in EPSG:3857. I've included a screenshot from QGIS showing OSM tiles, with the Lambert2008 geotiff on top. The pink contours are generated by pyhgtmap from the Lambert2008 tiff file, with the brown ones generated from the EPSG:3857 geotiff file.

Here's the listgeo output of that Lambert2008 geotiff file:

   Version: 1
   Key_Revision: 1.0
   Tagged_Information:
      ModelTiepointTag (2,3):
         0                 0                 0                
         323891.028851897  1030258.88565706  0                
      ModelPixelScaleTag (1,3):
         90                90                0                
      End_Of_Tags.
   Keyed_Information:
      GTModelTypeGeoKey (Short,1): ModelTypeProjected
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      GTCitationGeoKey (Ascii,8): "unknown"
      GeographicTypeGeoKey (Short,1): User-Defined
      GeogCitationGeoKey (Ascii,111): "GCS Name = unknown|Datum = Unknown based on GRS 1980 ellipsoid using towgs84=0,0,0,0,0,0,0|Primem = Greenwich|"
      GeogGeodeticDatumGeoKey (Short,1): User-Defined
      GeogAngularUnitsGeoKey (Short,1): Angular_Degree
      GeogEllipsoidGeoKey (Short,1): Ellipse_GRS_1980
      GeogSemiMajorAxisGeoKey (Double,1): 6378137          
      GeogInvFlatteningGeoKey (Double,1): 298.257222101    
      GeogPrimeMeridianLongGeoKey (Double,1): 0                
      GeogTOWGS84GeoKey (Double,3): 0                0                0                
      ProjectedCSTypeGeoKey (Short,1): User-Defined
      ProjectionGeoKey (Short,1): User-Defined
      ProjCoordTransGeoKey (Short,1): CT_LambertConfConic_2SP
      ProjLinearUnitsGeoKey (Short,1): Linear_Meter
      ProjStdParallel1GeoKey (Double,1): 49.8333333333333 
      ProjStdParallel2GeoKey (Double,1): 51.1666666666667 
      ProjFalseOriginLongGeoKey (Double,1): 4.35921583333333 
      ProjFalseOriginLatGeoKey (Double,1): 50.797815        
      ProjFalseOriginEastingGeoKey (Double,1): 649328           
      ProjFalseOriginNorthingGeoKey (Double,1): 665262           
      End_Of_Keys.
   End_Of_Geotiff.

Projection Method: CT_LambertConfConic_2SP
   ProjFalseOriginLatGeoKey: 50.797815 ( 50d47'52.13"N)
   ProjFalseOriginLongGeoKey: 4.359216 (  4d21'33.18"E)
   ProjStdParallel1GeoKey: 49.833333 ( 49d50' 0.00"N)
   ProjStdParallel2GeoKey: 51.166667 ( 51d10' 0.00"N)
   ProjFalseEastingGeoKey: 649328.000000 m
   ProjFalseNorthingGeoKey: 665262.000000 m
Ellipsoid: 7019/GRS 1980 (6378137.00,6356752.31)
TOWGS84: 0,0,0
Projection Linear Units: 9001/metre (1.000000m)

Corner Coordinates:
Upper Left    (  323891.029, 1030258.886)  (  0d35'44.32"W, 53d58'42.19"N)
Lower Left    (  323891.029,  353998.886)  (  0d 0'25.44"E, 47d54'50.38"N)
Upper Right   (  771911.029, 1030258.886)  (  6d13'40.63"E, 54d 3'43.06"N)
Lower Right   (  771911.029,  353998.886)  (  6d 0' 0.58"E, 47d59'15.17"N)
Center        (  547901.029,  692128.886)  (  2d54'47.76"E, 51d 1'49.64"N)

[agrenott]

Hi. Thanks for the bug report. Can you please attach/provide a link to the various tif files you used to reproduce the issue?

[rmanhaeve]

Here you go: lambert.tif 3857.tif

[agrenott]

The issue is due to the fact the projection used results in a "non rectangular" (there's probably a better denomination) area:

gdalinfo ~/Downloads/lambert.tif | grep Left 
Upper Left  (  323891.029, 1030258.886) (  0d35'44.32"W, 53d58'42.19"N)
Lower Left  (  323891.029,  353998.886) (  0d 0'25.44"E, 47d54'50.38"N)

Most of pyhgtmap internals make the strong assumption the area is rectangular, which I won't change. I'll add a check to detect such input and fail with an error instead.

[rmanhaeve]

Ok, thanks for looking into it!