Tuesday 28 January 2020

gdal - gdalwarp leaves horizontal artifacts regridding from EASE-Grid (laea) to Polarstero (stere)


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

Answer



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

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...