Invalid VRT Size when converting from EPSG:4326 to EPSG:3857

[vincentsarago]

ref: cogeotiff/rio-tiler#654

Expected behavior and actual behavior.

I've got a file which span the entire world (in epsg:4326, -180 90 180 -90), when building a VRT in EPSG:3857 the result size seems wrong.

I'm pretty sure this is because the bounds over flows the area of use of the EPSG:3857 projection.

gdalwarp -of VRT -t_srs epsg:3857 world_cog.tif world_cog.vrt
Creating output file that is 66P x 802L.

The Y size (802) seems ok but the X (66) is definitely wrong.

When I use dataset which doesn't cross the area of use everything is fine!

gdal_create -of COG -outsize 720 360 -a_srs EPSG:4326 -a_ullr -175 85 175 -85 world_cog.tif 
0...10...20...30...40...50...60...70...80...90...100 - done.

gdalwarp -of VRT -t_srs epsg:3857 world_cog.tif world_cog.vrt -overwrite
Creating output file that is 562P x 576L.
Processing world_cog.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

I'm not sure exactly if the issue in either in proj or in GDAL and if anything is possible 🤷

Steps to reproduce the problem.

gdal_create -of COG -outsize 720 360 -a_srs EPSG:4326 -a_ullr -180 90 180 -90 world_cog.tif 
0...10...20...30...40...50...60...70...80...90...100 - done.

gdalwarp -of VRT -t_srs epsg:3857 world_cog.tif world_cog.vrt
Creating output file that is 66P x 802L.
Processing world_cog.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

gdalinfo world_cog.vrt
Driver: VRT/Virtual Raster
Files: world_cog.vrt
       world_cog.tif
Size is 66, 802
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["unnamed",
        METHOD["Popular Visualisation Pseudo Mercator",
            ID["EPSG",1024]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (-20037508.342789243906736,242528680.943742722272873)
Pixel Size = (604620.393207252956927,-604620.393207252956927)
Metadata:
  AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  (-20037508.343,242528680.944) (180d 0' 0.00"W, 90d 0' 0.00"N)
Lower Left  (-20037508.343,-242376874.408) (180d 0' 0.00"W, 90d 0' 0.00"S)
Upper Right (19867437.609,242528680.944) (178d28'20.02"E, 90d 0' 0.00"N)
Lower Right (19867437.609,-242376874.408) (178d28'20.02"E, 90d 0' 0.00"S)
Center      (  -85035.367,   75903.268) (  0d45'49.99"W,  0d40'54.60"N)
Band 1 Block=66x128 Type=Byte, ColorInterp=Gray
  Overviews: 33x401

Operating system

Mac OS / CentOS (AWS Lambda)

GDAL version and provenance

GDAL 3.7.2 (homebrew)
GDAL 3.8 (docker)

gdalinfo --version
GDAL 3.8.0, released 2023/11/06

[jratike80]

Do you expect that gdalwarp would automatically crop the output according to the valid area of the target srs, in this case [World between 85.06°S and 85.06°N]? Like what can to done manually

gdalwarp -of VRT -t_srs epsg:3857 world_cog.tif world_cog.vrt -te -180 -85 180 85 -te_srs epsg:4326
Creating output file that is 699P x 697L.

Automatic clip feels like a heavy tool because sometimes it is OK to go over the area of use, for example if the source data comes around the boundary between adjacent UTM zones.

The problem comes from you asking gdalwarp to do something that is not possible, but how would you like to be informed about that?

[rouault]

I believe gdalwarp is behaving as it is designed. Latitude 90 projects to infinity in EPSG:3857 in theory. In practice due to numeric imprecision it converts to a huge northing value. Then gdalwarp tries to figure out an output resolution for that, which gives a huge skew toward the Y axis, and thus explains the contraction on the X axis
I don't see any potential "fix"

[vincentsarago]

Then gdalwarp tries to figure out an output resolution for that, which gives a huge skew toward the Y axis, and thus explains the contraction on the X axis

👍

I don't see any potential "fix"

That's what I thought too. Do you thing a warning will be possible for this case?

[jratike80]

It is certainly possible but is it worth it? Generally it is normal to go a little bit beyond the "area of use" and the transformation mathematics are accurate enough. So how to define the threshold for what is acceptable and what not? That said, you are not the first person that I have helped with this "whole world into Web Mercator" issue.

[vincentsarago]

It is certainly possible but is it worth it

Well right now this will lead to a really weird VRT

actual behaviour

expected behaviour

Just for context, I do not have an issue with a specific file but I'm just trying to make sure my library can handle this kind of edge cases.

At the end of the day, I do not really need a full VRT but I need the VRT resolution. In cogeotiff/rio-tiler#656 I'm try to get the resolution by reprojecting the bbox of the pixel at the center of the dataset... but I'm pretty sure there should be a better way for this.

[tbonfort]

I don't believe there's a clean solution to your problem here Vincent, getting the resolution requires dividing a geographical extent by some number of pixels, but in your case that extent is infinity 🤷‍♂️

[jratike80]

Right, it looks odd. It is caused by an user error but users do not know automatically what are the limits of the CRS what they are using. I am also just a user but I know the limits of 3857 because I have done that error so many times.

If EPSG:4326 -> EPSG:3857 is a common use case for your library, maybe you could implement the latitude cut -85 - 85 degrees into the library? But of course North Pole as EPSG:3857 is not the only impossible transformation that exists.

[rouault]

in #8740, I've added a heuristics for that special case. Hopefully I'm not opening a can of worms here...

[rcoup]

@rouault I was just writing a similar suggestion wrt EPSG:3857 since it's such a common target CRS.

I guess a more generic path would be Proj having a mechanism to warn users that they're operating outside the area of use for their CRS and they should expect odd things to happen.

[jratike80]

Actually that kind of mechanism is already running somewhere. I found 32 hits from gis.stackexchange about warnings like this

Warning 1: Unable to compute source region for output window 0,0,5789,4748, skipping.
ERROR 1: Too many points (439 out of 441) failed to transform,

Problem with a straight-lined Proj mechanism would be to decide the tolerances. Working outside the official area of use is rather common (e.g. combining data from adjacent UTM zones) and sometimes the transformations give reasonable and reversible results far away from the nominal area of use.

[rouault]

I guess a more generic path would be Proj having a mechanism to warn users that they're operating outside the area of use for their CRS and they should expect odd things to happen.

That's an interesting idea, but as often with map projections, things are tricky. For example for UTM projections, the recommended area of use is +/- 3 deg of longitude around the central meridian. But the ellipsoidal formulation of transverse mercator exists for all points of the spheroid, and the implementation in PROJ >= 4.9 is quite exact. Not fully exact, but it gives accurate results up to +/- 120 degree (perhaps even 150) away from the central meridian, and doing "gdalwarp autotest/gdrivers/data/small_world.tif out.tif -t_srs EPSG:32631 -overwrite" gives a very reasonable output for the whole world. So for a UTM projection, the area of use is mostly a loose indication, whereas for a Mercator one, it is close to be a constraint.

[rcoup]

That's an interesting idea, but as often with map projections, things are tricky.

Of course 😁 And the EPSG database isn't perfect either, I've personally filed a number of area-of-use bugfixes. I guess I was suggesting more of a method of warning that GDAL/OGR could propagate to the user.

$ gdalwarp -of VRT -t_srs epsg:3857 world_cog.tif world_cog.vrt
WARNING: coordinate transform is operating outside the nominal area of use for EPSG:3857, which may cause problems. See https://proj.org/E12345 for more details.
Creating output file that is 66P x 802L.

As you saying, ignoring the 6° bands for Transverse Mercator projections won't do anything violently terrible, but it might be a useful prompt for the user. We definitely couldn't do any automated clipping for most CRS (Polar Stereographic comes to mind alongside Mercator).