[GRASS-user] Re: NetCDF rotated pole

[Hermann Peifer]

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


=========================================================
How to re-project "rotated pole lat/lon grids" with GRASS
=========================================================

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.