Hi,
Is there any way to read a NetCDF file which has a rotated pole Mercator projection and have it inside GRASS with the North pole at North?
I'm trying to use AVISO SSH products on GRASS 6.4.1 and gdal 1.8.0
Thanks!
--
Eduardo Klein
Lab Sensores Remotos e
Instituto de Tecnología y Ciencias Marinas
Universidad Simón Bolívar
Caracas, Venezuela
ph (58) (212) 906-3111 ext 6700
fax (58) (212( 906-3111 ext 6701
On 19/03/2012 18:50, Eduardo Klein wrote:
Hi,
Is there any way to read a NetCDF file which has a rotated pole Mercator
projection and have it inside GRASS with the North pole at North?I'm trying to use AVISO SSH products on GRASS 6.4.1 and gdal 1.8.0
Thanks!
This is possible but works in a somewhat counter-intuitive way. See my below explanations which I have sent to a colleague, a couple of weeks ago. I assume that you can apply the steps in analogy.
Hope this helps, Hermann
0. On the command line, one can re-project from "rotated lat/lon coordinates" to "normal lat/lon coordinates" by using:
$ proj -v -f "%.6f" -m 57.2957795130823 +proj=ob_tran +o_proj=latlon +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180
...and for the opposite direction (from normal to rotated coordinates), one uses:
$ invproj -v -f "%.6f" -m 57.2957795130823 +proj=ob_tran +o_proj=latlon +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180
I noted that the proj/invproj (i.e. forward/inverse projection) logic is opposite to what I would expect.
1. gdalwarp cannot be used for re-projection, as GDAL/OGR only supports a subset of the projections known in PROJ.4. ob_tran is NOT part of the supported subset :-(.
2. GRASS *can* be used for re-projection, by doing the following steps:
a) Create a regular WGS84 location in GRASS, which has these projection parameters
==> PROJ_INFO <==
name: Latitude-Longitude
datum: wgs84
towgs84: 0.000,0.000,0.000
proj: ll
ellps: wgs84
==> PROJ_UNITS <==
unit: degree
units: degrees
meters: 1.0
b) Import the NetCDF data into WGS84 GRASS location, disable coordinate system check with -o
c) Create an ob_tran GRASS location, which has these projection parameters:
==> PROJ_INFO <==
name: General Oblique Transformation
datum: wgs84
towgs84: 0.000,0.000,0.000
proj: ob_tran
o_proj: latlon
ellps: wgs84
a: 6378137.0000000000
es: 0.0066943800
f: 298.2572235630
lat_0: 0.0000000000
lon_0: 180.0000000000
o_lat_p: 39.2500000000
o_lon_p: -162.0000000000
==> PROJ_UNITS <==
unit: meter
units: meters
meters: 0.0174532925
(The conversion factor 0. 0174532925 is pi/180, used for the degree/radians conversion.)
Example re-projection, from wgs84 to ob_tran location:
$ r.proj precip.1 location=wgs84 mapset=PERMANENT out=precip.1 method=nearest # or bilinear
NB: The procedure intentionally ignores that the original (rotated) lat/lon coordinates are based on a spherical Earth and just uses the given coordinate values on the WGS84 ellipsoid. This is usually better than trying to apply a "datum shift" between sphere and ellipsoid. The PROJ software behaves by default in a similar way.