5

I'm trying to understand how to use gdalwarp to warp an image. I believe I'm doing the basics correctly, but I might be missing some gdalwarp options.

The basic problem is that I see horizontal artifiacts in my output image.

Here's basic steps to reproduce.

Start with a simple 722x722 image. Here's one on imagur (Direct link to png)

Now use gdal_translate to apply metadata to make this a geotiff (this is a EASE-grid Lambert Azimuthal Equal Area projection on a 1924 authallic sphere.)

gdal_translate -a_srs '+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +a=6371228 +b=6371228 +units=m +no_defs' \
    -a_ullr -4524688.262500000 4524688.262500000 4524688.262500000 -4524688.262500000  \
    -co COMPRESS=DEFLATE -co PREDICTOR=2 -co ZLEVEL=9 \
    -of GTiff ./UiMbqSd.png ./UiMbqSd.withmetadata.tif

This seems to have the expected geographic information (~12.5km grid) confirmed with gdalinfo

 gdalinfo UiMBqSd.withmetadata.tif


Driver: GTiff/GeoTIFF
Files: UiMbqSd.withmetadata.tif
Size is 722, 722
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["unnamed ellipse",
        DATUM["unknown",
            SPHEROID["unnamed",6371228,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Azimuthal_Equal_Area"],
    PARAMETER["latitude_of_center",90],
    PARAMETER["longitude_of_center",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-4524688.262500000186265,4524688.262500000186265)
Pixel Size = (12533.762500000000728,-12533.762500000000728)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-4524688.263, 4524688.263) (135d 0' 0.00"W, 29d42'45.71"N)
Lower Left  (-4524688.263,-4524688.263) ( 45d 0' 0.00"W, 29d42'45.71"N)
Upper Right ( 4524688.263, 4524688.263) (135d 0' 0.00"E, 29d42'45.71"N)
Lower Right ( 4524688.263,-4524688.263) ( 45d 0' 0.00"E, 29d42'45.71"N)
Center      (   0.0000000,   0.0000000) (180d 0' 0.00"E, 90d 0' 0.00"N)
Band 1 Block=722x2 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA
Band 2 Block=722x2 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA
Band 3 Block=722x2 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA
Band 4 Block=722x2 Type=Byte, ColorInterp=Alpha

Next I regrid it to a polarstero projection.

gdalwarp -t_srs "+proj=stere +lat_0=90 +lon_0=-45 +lat_ts=70 +ellps=WGS84 +datum=WGS84 +units=m" \
    ./UiMbqSd.withmetadata.tif ./regridded_lon0_-45.tif

The problem I'm seeing is horizontal artifacts that appear to protrude around the middle of the image. Here's the png representation of the problem (http://imgur.com/Nrx4ZoS)

I thought it might be something weird about my corners, but I see these artifacts with different lat_0 values.

gdalwarp -t_srs "+proj=stere +lat_0=90 +lon_0=0 +lat_ts=70 +ellps=WGS84 +datum=WGS84 +units=m" \
    ./UiMbqSd.withmetadata.tif ./regridded_lon0_0.tif

Here's an example with no rotation (http://imgur.com/uiiF9Ir)

I'm currently running this on a mac:

> gdalwarp --version
GDAL 1.11.1, released 2014/09/24

But I've tested and seen the same behavior on UBUNTU 12.04:

GDAL 1.10.1, released 2013/08/26
  • Great examples and test image. It was simple to repeat the test. I got similar results with GDAL 2.0-dev on Windows. If you won't get an answer from here you can try gdal-dev mailing list. – user30184 Oct 7 '14 at 18:16
4

Try adding the -et (error threshold) option with lower thresholds than the default (0.125). When I use "-et 0.01", the horizontal artifacts disappear:

gdalwarp -t_srs "+proj=stere +lat_0=90 +lon_0=-45 +lat_ts=70 +ellps=WGS84 +datum=WGS84 +units=m" \
    -et 0.01 \
    ./UiMbqSd.withmetadata.tif ./regridded_lon0_-45.tif
  • I'm using gdal_translate and have the same issue, but there seems to be on -et option for gdal_translate. Any suggestions on what to do there? Thanks! – mike.bronner Mar 13 at 20:55
0

Try changing the resampling method for gdalwarp, e.g.

gdalwarp -r mode etc. etc. 

which seems to remove your problem. What makes most sense will be different for different data.

  • I just tried this with: gdalwarp -r {TYPE} -t_srs "+proj=stere +lat_0=90 +lon_0=0 +lat_ts=70 +ellps=WGS84 +datum=WGS84 +units=m" ./UiMbqSd.withmetadata.tif ./regridded_{TYPE}.tif Where TYPE = [average, lanczos, cubicspline, cubic, mode, bilinear, near] All of them continue to show the artifacts. For the record, the actual data I'm trying to regrid is ice age data, so it's coded integers and I need to keep the values from the input grid. Thanks Matt – Matt Savoie Oct 7 '14 at 21:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.