2

I've been working with the altimetry data from topex.ucsd.edu/WWW_html/mar_grav.html

Was having difficulties extracting the data, so I wrote a simple tool to do that. Creates a bitmap form the data source (result http://imgur.com/MV3ivxl)

The problem is: I want an image with rounded corners (is there a name for that?) rounded corners

but when I use gdalwarp, I'm getting a rectangle rectangle.

(actual result from gdalwarp command below: http://imgur.com/PZ2TEOe so you can see the projection is changing)

Trying to go from Mercator to Robinson:

# data readme says source is spherical Mercator, from 0,80.738 to 360,-80.738
$ gdal_translate -a_srs '+proj=merc +datum=WGS84' -a_ullr 0 80.738 360 -80.738 grav.bmp grav.tif
# save http://spatialreference.org/ref/esri/54030/ogcwkt/ to robinson_ogcwkt.txt
$ gdalwarp --debug on -wo SOURCE_EXTRA=125 -t_srs robinson_ogcwkt.txt grav.tif grav_robinson.tif

output from last command

GDAL_netCDF: 
=====
Open(), filename=[grav_robinson.tif]
GDAL_netCDF: 
=====
Open(), filename=[grav_robinson.tif]
GDAL: GDALOpen(grav.tif, this=0x1e00740) succeeds as GTiff.
OGRCT: PROJ >= 4.8.0 features enabled
OGRCT: Source: +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
OGRCT: Target: +proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
OGRCT: Source: +proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
OGRCT: Target: +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
OGRCT: Source: +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
OGRCT: Target: +proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
OGRCT: Source: +proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
OGRCT: Target: +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
Creating output file that is 6157P x 3148L.
GDAL: GDALDriver::Create(GTiff,grav_robinson.tif,6157,3148,3,Byte,(nil))
Processing input file grav.tif.
WARP: Copying metadata from first source to destination dataset
GDAL: GDALDefaultOverviews::OverviewScan()
GDAL: GDALWarpKernel()::GWKNearestNoMasksByte()
Src=0,0,2825x4320 Dst=0,0,3078x3148
0...10...20...30...40...GDAL: Potential thrashing on band 1 of grav.tif.
GDAL: GDALWarpKernel()::GWKNearestNoMasksByte()
Src=2574,0,2826x4320 Dst=3078,0,3079x3148
50...60...70...80...90...100 - done.
GDAL: Potential thrashing on band 1 of grav_robinson.tif.
GDAL: GDALClose(grav.tif, this=0x1e00740)
GDAL: 8640 block reads on 4320 block band 1 of grav.tif.
GDAL: GDALClose(grav_robinson.tif, this=0x1e05610)
GDAL: 12592 block reads on 3148 block band 1 of grav_robinson.tif.

additional info:

$ gdalinfo grav.tif 
Driver: GTiff/GeoTIFF
Files: grav.tif
Size is 5400, 4320
Coordinate System is:
PROJCS["unnamed",
    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"]],
    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"]]]
Origin = (0.000000000000000,80.738000000000000)
Pixel Size = (0.066666666666667,-0.037378703703704)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (   0.0000000,  80.7380000) (  0d 0' 0.01"E,  0d 0' 2.63"N)
Lower Left  (   0.0000000, -80.7380000) (  0d 0' 0.01"E,  0d 0' 2.63"S)
Upper Right (     360.000,      80.738) (  0d 0'11.64"E,  0d 0' 2.63"N)
Lower Right (     360.000,     -80.738) (  0d 0'11.64"E,  0d 0' 2.63"S)
Center      ( 180.0000000,   0.0000000) (  0d 0' 5.82"E,  0d 0' 0.01"N)
Band 1 Block=5400x1 Type=Byte, ColorInterp=Red
Band 2 Block=5400x1 Type=Byte, ColorInterp=Green
Band 3 Block=5400x1 Type=Byte, ColorInterp=Blue

$ gdalinfo grav_robinson.tif 
Driver: GTiff/GeoTIFF
Files: grav_robinson.tif
Size is 6157, 3148
Coordinate System is:
PROJCS["World_Robinson",
    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"]],
    PROJECTION["Robinson"],
    PARAMETER["longitude_of_center",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
    AUTHORITY["EPSG","9001"]]]
Origin = (0.000000000000000,78.093065584271017)
Pixel Size = (0.049619632572081,-0.049619632572081)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (   0.0000000,  78.0930656) (  0d 0' 0.01"E,  0d 0' 2.63"N)
Lower Left  (   0.0000000, -78.1095378) (  0d 0' 0.01"E,  0d 0' 2.63"S)
Upper Right (     305.508,      78.093) (  0d 0'11.64"E,  0d 0' 2.63"N)
Lower Right (     305.508,     -78.110) (  0d 0'11.64"E,  0d 0' 2.63"S)
Center      ( 152.7540389,  -0.0082361) (  0d 0' 5.82"E,  0d 0' 0.00"S)
Band 1 Block=6157x1 Type=Byte, ColorInterp=Red
Band 2 Block=6157x1 Type=Byte, ColorInterp=Green
Band 3 Block=6157x1 Type=Byte, ColorInterp=Blue
1

The error you are doing is that -a_ullr is in degrees, and not in the merc projection you are applying in the same step.

This will work:

gdal_translate -a_srs merc.txt -a_ullr -20037508 15987768 20037508 -15987768 grav.jpg grav.tif
gdalwarp --debug on -wo SOURCE_EXTRA=125 -t_srs robinson_ogcwkt.txt grav.tif grav_robinson.tif

enter image description here

Another thing you have to care about is that the bitmap ranges from 0° to 360°, while the robinson projection is from -180° to +180°. So the coordinates you get are actually shifted by 180°.

2

another tool for doing this is this superb freeware: http://www.giss.nasa.gov/tools/gprojector/ (for win, mac ana linux): more than 80 projections, just a few simple clicks, do yourself a favour and give it a try

  • Sadly that link is down – John Powell May 13 '14 at 19:44
  • @John Barca please try again - its online just now – Kurt May 13 '14 at 20:31
  • Wow, a lot of projections there, and I consider myself a projections geek. – John Powell May 13 '14 at 20:38

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.