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

Joaquim Luis jmfluis at gmail.com
Wed May 8 11:01:51 PDT 2019


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


More information about the PROJ mailing list