[GIS] gdalwarp results in incorrect extent

gdalgdalwarpgeotiff-tiff

I have a GeoTIFF that's 172800×86400 on the MODIS "250 meter" Sinusoidal grid (more like 231.7, but who's counting?). (You can find its metadata at the bottom of this post.) This bounty of resolution has turned out to be quite a pain. My goal is to get it into a global half-degree grid (720×360) for compatibility with other datasets I'm using.

The best solution appears to be the following. Steps 2 and 3 are based on info gleaned from this thread:

  1. Use gdalwarp to re-project into WGS84 using nearest-neighbor, with size 172800×86400, longitude range [-180 180], and latitude range [-90 90]. This will result in a map at 7.5-second resolution.

  2. Use gdaladdo to create a "240" overview using "average." This will result in a 720×360 image at 0.5-degree resolution.

  3. Use gdaltranslate to turn that overview into a GeoTIFF with all the proper tags and such.

Sounds good, right? (Of course if there's a better way, please let me know!) But when I run gdalwarp like so (without \ line breaks, which I added here for readability):

gdalwarp -s_srs "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs" \
-t_srs "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0" \
-te −180 −90 180 90 -ts 172800 86400 -co COMPRESS=PACKBITS infile.tif  outfile.tif

The result is a GeoTIFF (metadata at end of post) with the correct size, but the wrong corner coordinates. Instead of the entire world, it only contains the Northern part of the Eastern hemisphere (i.e., the quadrant of the Earth with positive longitudes and latitudes). The resolution is, correspondingly, ¼ of what it should be (0.25 degrees per side instead of 0.5).

I have tried various other structures of the gdalwarp command. These include:

  • Not specifying s_srs, which I got from the original's metadata (same result)
  • Specifying t_srs using the EPSG code (same result)
  • Combination of the above (same result)
  • Specifying -tr 0.002083333333333 0.002083333333333 -tap instead of -te and -ts. (I'm not sure what precision actually gets accepted into gdalwarp; that might be handy to know, too.)
    • This results in a GeoTIFF of the right resolution, but wrong latitudinal size/extent (172800×86350, slightly > 3 minutes short of the poles). One solution would be to somehow "pad" the rest of this with zeros, which is realistic for my data, to make it 172800×86400.

I know the other two steps work, as I went ahead and did them before realizing what gdalwarp had done. Like I said, let me know if there's a better way to do this, but I would really like to understand why I'm getting these results from gdalwarp.

In case it helps, I'm running Mac OS X 10.7.4 with GDAL 1.8.1. Thanks in advance!

Here is the metadata from gdalinfo for the original GeoTIFF (infile.tif):

Size is 172800, 86400
Coordinate System is:
PROJCS["SIN         E019",
    GEOGCS[,
        DATUM["unknown",
            SPHEROID["unretrievable - using WGS84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Sinusoidal"],
    PARAMETER["longitude_of_center",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-20015108.640000000596046,10007554.320000000298023)
Pixel Size = (231.656350000000003,-231.656350000000003)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=PACKBITS
  INTERLEAVE=BAND
Corner Coordinates:
ERROR 1: tolerance condition error
Upper Left  (-20015108.640,10007554.320)
ERROR 1: tolerance condition error
Lower Left  (-20015108.640,-10007554.320)
ERROR 1: tolerance condition error
Upper Right (20015108.640,10007554.320)
ERROR 1: tolerance condition error
Lower Right (20015108.640,-10007554.320)
Center      (   0.0000000,   0.0000000) (  0d 0' 0.01"E,  0d 0' 0.01"N)
Band 1 Block=172800x1 Type=Byte, ColorInterp=Gray

And here's the metadata from the result of gdalwarp (outfile.tif):

Size is 172800, 86400
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 = (0.000000000000000,90.000000000000000)
Pixel Size = (0.001041666666667,-0.001041666666667)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=PACKBITS
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   0.0000000,  90.0000000) (  0d 0' 0.01"E, 90d 0' 0.00"N)
Lower Left  (   0.0000000,   0.0000000) (  0d 0' 0.01"E,  0d 0' 0.01"N)
Upper Right ( 180.0000000,  90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"N)
Lower Right ( 180.0000000,   0.0000000) (180d 0' 0.00"E,  0d 0' 0.01"N)
Center      (  90.0000000,  45.0000000) ( 90d 0' 0.00"E, 45d 0' 0.00"N)
Band 1 Block=172800x1 Type=Byte, ColorInterp=Gray

Best Answer

I figured out the answer. Turns out the note-taking program I was using to compose the code before pasting into Terminal auto-corrected the hyphens/dashes in front of the first two -te elements. They got changed from simple hyphens (-180 -90) to en dashes (–180 –90). Changing those back seems to have fixed everything.

Related Question