使用gdalwarp从wgs84转换到EPSG:3857不覆盖世界

Question

我使用的是gdal 1.10和2.1.1。

我在WGS84中有一个VRT数据源，其中我将角坐标强制为EPSG的最小/最大值:3857(-180,85.5180，-85.5)。

此VRT的gdalinfo输出如下所示:

Size is 1296001, 601200
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-180.000000000000000,85.500000000000000)
Pixel Size = (0.000277777563443,-0.000284431137725)
Corner Coordinates:
Upper Left  (-180.0000000,  85.5000000) (180d 0' 0.00"W, 85d30' 0.00"N)
Lower Left  (-180.0000000, -85.5000000) (180d 0' 0.00"W, 85d30' 0.00"S)
Upper Right ( 180.0000000,  85.5000000) (180d 0' 0.00"E, 85d30' 0.00"N)
Lower Right ( 180.0000000, -85.5000000) (180d 0' 0.00"E, 85d30' 0.00"S)
Center      (   0.0000000,  -0.0000000) (  0d 0' 0.01"E,  0d 0' 0.00"S)
Band 1 Block=128x128 Type=Int16, ColorInterp=Gray

基本上，我有世界减去极点。

现在我想把它转换成EPSG:3857。

我使用gdalwarp，使用双线性插值:

./gdalwarp -of VRT -co TILED=YES -srcnodata 9999 -t_srs 'EPSG:3785'  -multi  wgs84.vrt  wmerc.vrt  -overwrite -r bilinear

在wmerc上运行gdalinfo会得到这个:

Size is 995026, 1025175
Coordinate System is:
PROJCS["Popular Visualisation CRS / Mercator (deprecated)",
    GEOGCS["Popular Visualisation CRS",
        DATUM["Popular_Visualisation_Datum",
            SPHEROID["Popular Visualisation Sphere",6378137,0,
                AUTHORITY["EPSG","7059"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6055"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4055"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],
    AUTHORITY["EPSG","3785"]]
Origin = (-20037508.342789247632027,20644642.363762538880110)
Pixel Size = (40.275313937642778,-40.275354414328227)
Corner Coordinates:
Upper Left  (-20037508.343,20644642.364) (180d 0' 0.00"E, 85d28'11.27"N)
Lower Left  (-20037508.343,-20644644.098) (180d 0' 0.00"E, 85d28'11.28"S)
Upper Right (20037476.183,20644642.364) (179d59'58.96"E, 85d28'11.27"N)
Lower Right (20037476.183,-20644644.098) (179d59'58.96"E, 85d28'11.28"S)
Center      ( -16.0797308,  -0.8670919) (  0d 0' 0.52"W,  0d 0' 0.03"S)
Band 1 Block=512x128 Type=Int16, ColorInterp=Gray
  NoData Value=9999

注意左上/左下的角坐标看起来是正确的，但是右上/右下的角坐标(经度)缺失了32个单位。广义地说，我在右边缺少一条银条，但只是在坐标方面。数据在那里，但右边的坐标似乎不对。

为什么?

i 可以只是修改坐标来匹配世界(纵向)使用gdal_translate，但我担心我在这里忽略了一些可能会回来咬我的东西。

Answer 1

EPSG 3857限定纬度在-85到80度之间。它这样做有两个原因:

1. 它允许整个地球是方形的，并通过四叉树进行索引，其中每个节点都是方形的。

2. 墨卡托投影在两极爆炸;考虑到大多数使用地图的人对地图不感兴趣，投影一举两得，抛弃了它们。