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

Wed May 8 09:44:32 PDT 2019

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 <schwehr at gmail.com> 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 <jmfluis at gmail.com> 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 <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
>>
>>
>>
>
>
> --
> --
> http://schwehr.org
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20190508/8909f443/attachment-0001.html>