[PROJ] Projection for Sandwell et al.'s topex topo and grav files?

Tue May 7 18:52:11 PDT 2019

The IMG grids are spherical Mercator inches so there is likely no EPSG or WKT that represents that. But you can do all the work in GMT.

1- Convert to netcdf, maintaining the Merc projection, but change the origin to (0,0) (it was in the LL corner). That’s the role of -C

img2grd -R-180/180/-80/80 topo_19.1.img -Gspherical-mercator-proj.grd -T1 -D -M -C

2- Edit the header to convert the inches to Mercator meters. The grid above has x_min = -180; x_max = 180.

And that corresponds to 6378137 *2pi = 4.007501668557849e7  meters. So we can compute a scale factor of 

111319.49079327358 = 6378137 *2pi / 360

and apply it to the region in degrees.

[-180 180 -139.6 139.6] .* 111319.49079327358 = -20037508.3427892 20037508.3427892 -15540200.914741 15540200.914741

So finaly use grdedit to change the limits and assign it a proj4 string describing the projection

grdedit spherical-mercator-proj.grd -R-20037508.3427892/20037508.3427892/-15540200.914741/15540200.914741 -J"+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m +no_defs"

You can now confirm with gdalinfo

C:\v>gdalinfo spherical-mercator-proj.grd

Warning 1: dimension #1 (x) is not a Longitude/X dimension.

Warning 1: dimension #0 (y) is not a Latitude/Y dimension.

Driver: netCDF/Network Common Data Format

Files: spherical-mercator-proj.grd

Size is 21600, 16752

Coordinate System is:

PROJCS["unnamed",

    GEOGCS["unnamed ellipse",

        DATUM["unknown",

            SPHEROID["unnamed",6378137,0]],

        PRIMEM["Greenwich",0],

        UNIT["degree",0.0174532925199433]],

    PROJECTION["Mercator_1SP"],

    PARAMETER["central_meridian",0],

    PARAMETER["scale_factor",1],

    PARAMETER["false_easting",0],

    PARAMETER["false_northing",0],

    UNIT["Meter",1]]

Origin = (-20037148.342789199000000,15540200.914741000000000)

Pixel Size = (1855.324846554555500,-1855.324846554560700)

Metadata:

  grid_mapping#spatial_ref=PROJCS["unnamed",

    GEOGCS["unnamed ellipse",

        DATUM["unknown",

            SPHEROID["unnamed",6378137,0]],

        PRIMEM["Greenwich",0],

        UNIT["degree",0.0174532925199433]],

    PROJECTION["Mercator_1SP"],

    PARAMETER["central_meridian",0],

    PARAMETER["scale_factor",1],

    PARAMETER["false_easting",0],

    PARAMETER["false_northing",0],

    UNIT["Meter",1]]

  NC_GLOBAL#Conventions=CF-1.7

  NC_GLOBAL#description=Spherical Mercator Projected with -Jm1 -R-180/180/-80.0023237126/80.0023237126

  NC_GLOBAL#GMT_version=6.0.0_bcb87fa-dirty_2019.05.07 [64-bit] [MP]

  NC_GLOBAL#history=img2grd -R-180/180/-80/80 topo_19.1.img -Gspherical-mercator-proj.grd -T1 -D -M -C

  NC_GLOBAL#node_offset=1

  NC_GLOBAL#title=Data from Altimetry

  x#actual_range={-20037148.3427892,20037868.3427892}

  x#long_name=x_units

  y#actual_range={-15540200.914741,15540200.914741}

  y#long_name=y_units

  z#actual_range={-10914,8550}

  z#grid_mapping=grid_mapping

  z#long_name=meters, mGal, Eotvos, micro-radians or Myr, depending on img file and -S.

  z#_FillValue=-1.#IND

Corner Coordinates:

Upper Left  (-20037148.343,15540200.915) (179d59'48.36"W, 80d 0' 8.37"N)

Lower Left  (-20037148.343,-15540200.915) (179d59'48.36"W, 80d 0' 8.37"S)

Upper Right (20037868.343,15540200.915) (179d59'48.36"W, 80d 0' 8.37"N)

Lower Right (20037868.343,-15540200.915) (179d59'48.36"W, 80d 0' 8.37"S)

Center      (     360.000,       0.000) (  0d 0'11.64"E,  0d 0' 0.01"N)

Band 1 Block=21600x1 Type=Float32, ColorInterp=Undefined

  NoData Value=nan

  Metadata:

    actual_range={-10914,8550}

    grid_mapping=grid_mapping

    long_name=meters, mGal, Eotvos, micro-radians or Myr, depending on img file and -S.

    NETCDF_VARNAME=z

    _FillValue=-1.#IND

Joaquim

(with Paul’s help for the scaling factor)

From: PROJ <proj-bounces at lists.osgeo.org> On Behalf Of Kurt Schwehr
Sent: Tuesday, May 7, 2019 10:52 PM
To: PROJ <proj at lists.osgeo.org>
Subject: [PROJ] Projection for Sandwell et al.'s topex topo and grav files?

Hi all,

I figured I should ask here too if anyone know what the correct projection is for the Sandwell .img grids from http://topex.ucsd.edu.  I'm trying to keep the files in their original projections as I switch them to geotiffs.  I've asked David Sandwell directly too if he knows.

This works, but warps the data first:

gmt img2grd -V -R-180/180/-80/80 topo_18.1.img  -Gtest2.grd -T1 -D

gdal_translate -a_srs EPSG:4326 test2.grd test2.tif

gdalinfo test2.tif  # Results look believable

gdal_translate topo-18-1-epsg4326.tif topo-18-1-epsg4326-deflate.tif -co COMPRESS=DEFLATE -co PREDICTOR=3

Then imported into QGIS or Earth Engine as a normal user, things line up pretty well.

    https://code.earthengine.google.com/f43c8b13bafa26fd8c7f83ce9a919f6e

I'd rather do it more like this:

gmt img2grd -V -R-180/180/-80/80 topo_18.1.img  -Gspherical-mercator-proj.grd -T1 -D -M

img2grd: Expects topo_18.1.img to be 21600 by 17280 pixels spanning 0/360.0/-80.738009/80.738009.

img2grd: To fit [averaged] input, your topo_18.1.img is adjusted to -R-180/180/-80.0023237126/80.0023237126.

img2grd: The output grid size will be 21600 by 16752 pixels.

img2grd: Created 21600 by 16752 Mercatorized grid file.  Min, Max values are -10914  8550

gives this which doesn't work as is:

gdalinfo spherical-mercator-proj.grd

Driver: netCDF/Network Common Data Format

Files: spherical-mercator-proj.grd

Size is 21600, 16752

Coordinate System is `'

Origin = (0.000000000000000,279.199999999999989)

Pixel Size = (0.016666666666667,-0.016666666666667)

Metadata:

  NC_GLOBAL#Conventions=COARDS, CF-1.5

  NC_GLOBAL#description=Spherical Mercator Projected with -Jm1 -R0/360/-80.0023237126/80.0023237126

  NC_GLOBAL#GMT_version=5.4.3 (r19528) [64-bit]

Should I be using one of 54004 or 41001?  e.g. https://epsg.io/?q=spherical+mercator

Thanks!

-kurt

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20190508/e6dbe2c6/attachment-0001.html>