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

classic Classic list List threaded Threaded
6 messages Options
Reply | Threaded
Open this post in threaded view
|

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

Kurt Schwehr-2
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.



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


_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Projection for Sandwell et al.'s topex topo and grav files?

J Luis

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 <[hidden email]> On Behalf Of Kurt Schwehr
Sent: Tuesday, May 7, 2019 10:52 PM
To: PROJ <[hidden email]>
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.

 

 

 

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

 


_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Projection for Sandwell et al.'s topex topo and grav files?

Kurt Schwehr-2
Thanks!  I'll be investigating more today.  Monica Schwehr mentioned that you have a different Earth radius than she was seeing in the GMT code base.  I'll follow up on that after I get a chance to look through GMT head and GMT 4 more.

On Tue, May 7, 2019 at 6:52 PM <[hidden email]> wrote:

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 <[hidden email]> On Behalf Of Kurt Schwehr
Sent: Tuesday, May 7, 2019 10:52 PM
To: PROJ <[hidden email]>
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.

 

 

 

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

 



--

_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Projection for Sandwell et al.'s topex topo and grav files?

J Luis

Right, I used a radius from a proj4 string that I had here but I'm realizing that 6378137 is WGS's 84 major axis. OK, I should have used the proj
 sphere a=6370997.0
(or the GMT's 6371008.7714, but that would make no difference as long as we calculate the scale factor with one and use that same one in the proj string)
With it, and recomputing the scale factor, you should use

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

However, I realized that grdedit is doing something to the limits that I need to understand better. See that it slightly changed the XX limits that we provided. x_min & x_max are no longer symetric, but y_min,max are.

C:\v>grdinfo spherical-mercator-proj.grd
spherical-mercator-proj.grd: Title: Data from Altimetry
spherical-mercator-proj.grd: Command: img2grd -R-180/180/-80/80 topo_19.1.img -Gspherical-mercator-proj.grd -T1 -D -M -C
spherical-mercator-proj.grd: Remark: Spherical Mercator Projected with -Jm1 -R-180/180/-80.0023237126/80.0023237126
spherical-mercator-proj.grd: Pixel node registration used [Cartesian grid]
spherical-mercator-proj.grd: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7
spherical-mercator-proj.grd: x_min: -20014717.3712 x_max: 20015437.3712 x_inc: 1853.24790474 name:  n_columns: 21600
spherical-mercator-proj.grd: y_min: -15522804.4501 y_max: 15522804.4501 y_inc: 1853.24790474 name:  n_rows: 16752


On Wed, May 8, 2019 at 4:11 PM Kurt Schwehr <[hidden email]> wrote:
Thanks!  I'll be investigating more today.  Monica Schwehr mentioned that you have a different Earth radius than she was seeing in the GMT code base.  I'll follow up on that after I get a chance to look through GMT head and GMT 4 more.

On Tue, May 7, 2019 at 6:52 PM <[hidden email]> wrote:

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 <[hidden email]> On Behalf Of Kurt Schwehr
Sent: Tuesday, May 7, 2019 10:52 PM
To: PROJ <[hidden email]>
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.

 

 

 

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

 



--

_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Projection for Sandwell et al.'s topex topo and grav files?

Kurt Schwehr-2
Thanks Joaquim.  

I'll see if I can get that in and visualized.  I'm having a challenge with alignment with a hill shaded background.

Life is extra fun that src/Ellipsoids.txt was replace by by gmt_ellipsoids.h back in January:

git log --all --full-history -- src/Ellipsoids.txt

I think I'm able to trace what is going on:

./src/gmt_ellipsoids.h: {"Sphere", 1984, 6371008.7714, 0},

egrep 'Jm|Sphere|grdproject' src/img/img2grd.c
 * preserving the Mercator projection (gmtdefaults PROJ_ELLIPSOID = Sphere)
 * and height of the (spherical) Mercator map with -Rww/ee/ss/nn and -Jm1.
/* Set up header with Mercatorized dimensions assuming -Jm1  */
snprintf (Merc->header->x_units, GMT_GRID_UNIT_LEN80, "Spherical Mercator projected Longitude, -Jm1, length from %.12g", left);
snprintf (Merc->header->y_units, GMT_GRID_UNIT_LEN80, "Spherical Mercator projected Latitude, -Jm1, length from %.12g", bottom);
snprintf (Merc->header->remark, GMT_GRID_REMARK_LEN160, "Spherical Mercator Projected with -Jm1 -R%.12g/%.12g/%.12g/%.12g", wesn[XLO], wesn[XHI], south2, north2);
GMT_Report (API, GMT_MSG_LONG_VERBOSE, "Undo the implicit spherical Mercator -Jm1i projection.\n");
/* Preparing source and destination for GMT_grdproject */
/* a. Register the Mercator grid to be the source read by GMT_grdproject by passing a pointer */
GMT_Report (API, GMT_MSG_DEBUG, "Open Mercator Grid as grdproject virtual input\n");
/* b. If -E: Register a grid struct Geo to be the destination allocated and written to by GMT_grdproject, else write to -G<file> */
GMT_Report (API, GMT_MSG_DEBUG, "Register memory Grid container as grdproject output\n");
sprintf (cmd, "-R%g/%g/%g/%g -Jm1i -I %s -G%s --PROJ_ELLIPSOID=Sphere --PROJ_LENGTH_UNIT=inch --GMT_HISTORY=false", west, east, south2, north2, input, output);
GMT_Report (API, GMT_MSG_DEBUG, "Calling grdproject %s.\n", cmd);
if (GMT_Call_Module (API, "grdproject", GMT_MODULE_CMD, cmd)!= GMT_NOERROR) { /* Inverse project the grid or fail */
/* a. Register the Geographic grid returned by GMT_grdproject to be the source read by GMT_grdsample by passing a pointer */
GMT_Report (API, GMT_MSG_DEBUG, "Read Geo Grid container as grdproject virtual output\n");

On Wed, May 8, 2019 at 9:44 AM J Luis <[hidden email]> wrote:

Right, I used a radius from a proj4 string that I had here but I'm realizing that 6378137 is WGS's 84 major axis. OK, I should have used the proj
 sphere a=6370997.0
(or the GMT's 6371008.7714, but that would make no difference as long as we calculate the scale factor with one and use that same one in the proj string)
With it, and recomputing the scale factor, you should use

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

However, I realized that grdedit is doing something to the limits that I need to understand better. See that it slightly changed the XX limits that we provided. x_min & x_max are no longer symetric, but y_min,max are.

C:\v>grdinfo spherical-mercator-proj.grd
spherical-mercator-proj.grd: Title: Data from Altimetry
spherical-mercator-proj.grd: Command: img2grd -R-180/180/-80/80 topo_19.1.img -Gspherical-mercator-proj.grd -T1 -D -M -C
spherical-mercator-proj.grd: Remark: Spherical Mercator Projected with -Jm1 -R-180/180/-80.0023237126/80.0023237126
spherical-mercator-proj.grd: Pixel node registration used [Cartesian grid]
spherical-mercator-proj.grd: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7
spherical-mercator-proj.grd: x_min: -20014717.3712 x_max: 20015437.3712 x_inc: 1853.24790474 name:  n_columns: 21600
spherical-mercator-proj.grd: y_min: -15522804.4501 y_max: 15522804.4501 y_inc: 1853.24790474 name:  n_rows: 16752


On Wed, May 8, 2019 at 4:11 PM Kurt Schwehr <[hidden email]> wrote:
Thanks!  I'll be investigating more today.  Monica Schwehr mentioned that you have a different Earth radius than she was seeing in the GMT code base.  I'll follow up on that after I get a chance to look through GMT head and GMT 4 more.

On Tue, May 7, 2019 at 6:52 PM <[hidden email]> wrote:

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 <[hidden email]> On Behalf Of Kurt Schwehr
Sent: Tuesday, May 7, 2019 10:52 PM
To: PROJ <[hidden email]>
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.

 

 

 

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

 



--


--

_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Projection for Sandwell et al.'s topex topo and grav files?

J Luis
Kurt, I’ll be away for a couple hours so cannot test this, but if the limits change in the grdedit step are wrong I suspect you can split into two steps:
1 change only the -R
2. Add the project info with the -J part. 

Sent from my iDedo

No dia 08/05/2019, às 18:50, Kurt Schwehr <[hidden email]> escreveu:

Thanks Joaquim.  

I'll see if I can get that in and visualized.  I'm having a challenge with alignment with a hill shaded background.

Life is extra fun that src/Ellipsoids.txt was replace by by gmt_ellipsoids.h back in January:

git log --all --full-history -- src/Ellipsoids.txt

I think I'm able to trace what is going on:

./src/gmt_ellipsoids.h: {"Sphere", 1984, 6371008.7714, 0},

egrep 'Jm|Sphere|grdproject' src/img/img2grd.c
 * preserving the Mercator projection (gmtdefaults PROJ_ELLIPSOID = Sphere)
 * and height of the (spherical) Mercator map with -Rww/ee/ss/nn and -Jm1.
/* Set up header with Mercatorized dimensions assuming -Jm1  */
snprintf (Merc->header->x_units, GMT_GRID_UNIT_LEN80, "Spherical Mercator projected Longitude, -Jm1, length from %.12g", left);
snprintf (Merc->header->y_units, GMT_GRID_UNIT_LEN80, "Spherical Mercator projected Latitude, -Jm1, length from %.12g", bottom);
snprintf (Merc->header->remark, GMT_GRID_REMARK_LEN160, "Spherical Mercator Projected with -Jm1 -R%.12g/%.12g/%.12g/%.12g", wesn[XLO], wesn[XHI], south2, north2);
GMT_Report (API, GMT_MSG_LONG_VERBOSE, "Undo the implicit spherical Mercator -Jm1i projection.\n");
/* Preparing source and destination for GMT_grdproject */
/* a. Register the Mercator grid to be the source read by GMT_grdproject by passing a pointer */
GMT_Report (API, GMT_MSG_DEBUG, "Open Mercator Grid as grdproject virtual input\n");
/* b. If -E: Register a grid struct Geo to be the destination allocated and written to by GMT_grdproject, else write to -G<file> */
GMT_Report (API, GMT_MSG_DEBUG, "Register memory Grid container as grdproject output\n");
sprintf (cmd, "-R%g/%g/%g/%g -Jm1i -I %s -G%s --PROJ_ELLIPSOID=Sphere --PROJ_LENGTH_UNIT=inch --GMT_HISTORY=false", west, east, south2, north2, input, output);
GMT_Report (API, GMT_MSG_DEBUG, "Calling grdproject %s.\n", cmd);
if (GMT_Call_Module (API, "grdproject", GMT_MODULE_CMD, cmd)!= GMT_NOERROR) { /* Inverse project the grid or fail */
/* a. Register the Geographic grid returned by GMT_grdproject to be the source read by GMT_grdsample by passing a pointer */
GMT_Report (API, GMT_MSG_DEBUG, "Read Geo Grid container as grdproject virtual output\n");

On Wed, May 8, 2019 at 9:44 AM J Luis <[hidden email]> wrote:

Right, I used a radius from a proj4 string that I had here but I'm realizing that 6378137 is WGS's 84 major axis. OK, I should have used the proj
 sphere a=6370997.0
(or the GMT's 6371008.7714, but that would make no difference as long as we calculate the scale factor with one and use that same one in the proj string)
With it, and recomputing the scale factor, you should use

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

However, I realized that grdedit is doing something to the limits that I need to understand better. See that it slightly changed the XX limits that we provided. x_min & x_max are no longer symetric, but y_min,max are.

C:\v>grdinfo spherical-mercator-proj.grd
spherical-mercator-proj.grd: Title: Data from Altimetry
spherical-mercator-proj.grd: Command: img2grd -R-180/180/-80/80 topo_19.1.img -Gspherical-mercator-proj.grd -T1 -D -M -C
spherical-mercator-proj.grd: Remark: Spherical Mercator Projected with -Jm1 -R-180/180/-80.0023237126/80.0023237126
spherical-mercator-proj.grd: Pixel node registration used [Cartesian grid]
spherical-mercator-proj.grd: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7
spherical-mercator-proj.grd: x_min: -20014717.3712 x_max: 20015437.3712 x_inc: 1853.24790474 name:  n_columns: 21600
spherical-mercator-proj.grd: y_min: -15522804.4501 y_max: 15522804.4501 y_inc: 1853.24790474 name:  n_rows: 16752


On Wed, May 8, 2019 at 4:11 PM Kurt Schwehr <[hidden email]> wrote:
Thanks!  I'll be investigating more today.  Monica Schwehr mentioned that you have a different Earth radius than she was seeing in the GMT code base.  I'll follow up on that after I get a chance to look through GMT head and GMT 4 more.

On Tue, May 7, 2019 at 6:52 PM <[hidden email]> wrote:

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 <[hidden email]> On Behalf Of Kurt Schwehr
Sent: Tuesday, May 7, 2019 10:52 PM
To: PROJ <[hidden email]>
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.

 

 

 

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

 



--


--

_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj