[PROJ] Projection for Sandwell et al.'s topex topo and grav files?
Kurt Schwehr
schwehr at gmail.com
Wed May 8 08:11:17 PDT 2019
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/cbc21120/attachment-0001.html>
More information about the PROJ
mailing list