Re: [CF-metadata] Rotated Pole projection possibly wrong

2013-03-09 Thread Tom Kunicki

The Z, Y', X'' rotation scheme, intrinsic or extrinsic, sounds like an 
application of Euler (sounds like Oiler) angles [1,2].

Tom Kunicki
Center for Integrated Data Analytics
U.S. Geological Survey
8505 Research Way
Middleton, WI  53562

[1] http://en.wikipedia.org/wiki/Euler_angles
[2] http://mathworld.wolfram.com/EulerAngles.html


On Mar 6, 2013, at 2:17 PM, Seth McGinnis mcgin...@ucar.edu wrote:

 So there are an infinite number of ways to decompose a 3-D orientation into
 three orthogonal rotations.  Specifying the projection in terms of the final
 location of the north pole and an optional rotation around that point means
 that CF is actually agnostic about exactly how you get there.
 
 That said, given that it's being specified in that way, the decomposition that
 maps most straightforwardly to those parameters is a Z-Y-X rotation sequence
 using an extrinsic (space-fixed) coordinate system:
 
 First rotate by angle phi around the north pole in space-fixed coordinates.
 Then rotate by angle theta around 0N 90E in space-fixed coordinates.
 Finally, rotate by angle psi around 0N 0E in space-fixed coordinates.
 
 The location of the grid's north pole in the original lat-lon coordinate 
 system
 is now at longitude = phi, colatitude = minus theta.  
 
 (I'm not sure how psi relates to north_pole_grid_longitude mathematically; it
 defaults to zero and I've never seen it in use, so I think it's included 
 mostly
 for completeness. I agree that we may well want to change it to something more
 obvious like an azimuth angle.)
 
 Switching to the GRIB2 convention where you use the south pole instead of the
 north should be just a matter of adding 180 degrees to phi and negating theta.
 
 Cheers,
 
 --Seth
 
 
 On Wed, 6 Mar 2013 18:29:47 +
 David Hassell d.c.hass...@reading.ac.uk wrote:
 Hello,
 
 It is a bit of mess! As I understand it, the full rotation described
 is a sequence of rotations about three different axes:
 
  Z  = 90S - 90N
  Y' = 90W' - 90E'
  X = 0E - 180E
 
 where it is understood that the definitions of the N-S and W-E axes
 change after each rotation (hence the primes and double
 primes). Therefore the order in which they are done matters.
 
 I suspect that the usual and assumed order is Z, Y', X?
 
 From the GRIB-2 stuff John posted the north_pole_grid_longitude gives
 the rotation about the X axis.
 
 Something like John's Angle of rotation seems right to me. 
 
 The apparent lack of consistency between the parameter names irks
 me. Perhaps one solution could be to:
 
 i) Add some text to the conventions state the order (Z, Y', X, say)
  and direction of rotations.
 
 ii) add three new, consistent, self-describing parameters (e.g.)
 
 angle_of_rotation_z
 angle_of_rotation_y
 angle_of_rotation_x
 
 iii) allow aliases for backwards compatibility
 
 grid_north_pole_longitude = angle_of_rotation_z
 grid_north_pole_latitude  = angle_of_rotation_y
 north_pole_grid_longitude = angle_of_rotation_x   
 
 Or would that just obfuscate things even more?
 
 All the best,
 
 David
 
 P.S. If you have a copy available, there are some nice descriptions in
Coordinate Systems and Map Projections by D. H. Maling
 
 --
 David Hassell
 National Centre for Atmospheric Science (NCAS)
 Department of Meteorology, University of Reading,
 Earley Gate, PO Box 243,
 Reading RG6 6BB, U.K.
 
 Tel   : 0118 3785613
 Fax   : 0118 3788316
 E-mail: d.c.hass...@reading.ac.uk
 ___
 CF-metadata mailing list
 CF-metadata@cgd.ucar.edu
 http://mailman.cgd.ucar.edu/mailman/listinfo/cf-metadata
 
 ___
 CF-metadata mailing list
 CF-metadata@cgd.ucar.edu
 http://mailman.cgd.ucar.edu/mailman/listinfo/cf-metadata

___
CF-metadata mailing list
CF-metadata@cgd.ucar.edu
http://mailman.cgd.ucar.edu/mailman/listinfo/cf-metadata


Re: [CF-metadata] Rotated Pole projection possibly wrong

2013-03-06 Thread John Caron

Hi Heiko:

So

grid_north_pole_longitude = 
normalizeLongitude180(longitudeOfSouthernPoleInDegrees)

grid_north_pole_latitude = -1 * longitudeOfSouthernPoleInDegrees

?

where does one find documentation on proj4's ob_tran routine?

thanks,
John



On 3/6/2013 8:03 AM, Heiko Klein wrote:

Hi John,

all our models use rotated_latitute_longitude. And we have a 
consistent way of translating it between grib, proj and netcdf-java, 
i.e. publicly available at:

http://thredds.met.no/thredds/catalog/metno/proff4km/default/catalog.html


grid_mapping_name: rotated_latitude_longitude
grid_north_pole_longitude: 156.0
grid_north_pole_latitude: 23.5
earth_radius: 6367470.0

to proj:
proj4: +proj=ob_tran +o_proj=longlat +lon_0=-24 +o_lat_p=23.5 
+a=6367470.0 +no_defs


from grib:
lonRot = longitudeOfSouthernPoleInDegrees
latRot = latitudeOfSouthernPoleInDegrees

+proj=ob_tran +o_proj=longlat +lon_0=  
normalizeLongitude180(lonRot)   +o_lat_p=  (-1 * latRot);




This is the set of parameters which have been useful. I'm not sure if 
they are named in a useful way. proj4's ob_tran even allows for a 
o_lon_p parameter, but this doesn't seem to be used by grib or CF.


Best regards,

Heiko



On 2013-03-06 15:18, John Caron wrote:

The Rotated Pole projection here:

http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.6/cf-conventions.html#appendix-grid-mappings 




Rotated pole

grid_mapping_name = rotated_latitude_longitude


Map parameters:

  *

grid_north_pole_latitude

  *

grid_north_pole_longitude

  *

north_pole_grid_longitude - This parameter is option (default 
is 0).


Map coordinates:

The rotated latitude and longitude coordinates are identified by
the standard_name attribute values
grid_latitude and grid_longitude respectively.



is problematic. It has no reference to Proj.4, and is not in Snyder.

Its likely the same as the GRIB-2 GDS Template 3.1:

Grid Definition Template 3.1: Rotated Latitude/longitude (or
equidistant cylindrical, or Plate Carrée)

Octet No. Contents

   15-72 Same as Grid Definition Template 3.0 (see Note 1)

   73-76 Latitude of the southern pole of projection

77-80 Longitude of the southern pole of
projection

   81-84 Angle of rotation of projection

   85-nnList of number of points along
each meridian or parallel (These octets are only present for
quasi-regular grids as described in Note 3)

Notes:


(2) Three parameters define a general latitude/longitude
coordinate system, formed by a general rotation of the sphere. One
choice for these parameters is:

(a) The geographic latitude in degrees of the southern pole of
the coordinate system, θp for example.

(b) The geographic longitude in degrees of the southern pole of
the coordinate system, λp for example.

(c) The angle of rotation in degrees about the new polar axis
(measured clockwise when looking from the southern to the northern pole)
of the coordinate system, assuming the new axis to have been obtained by
first rotating the sphere through λp degrees about the geographic polar
axis, and then rotating through (90 + θp) degrees so that the southern
pole moved along the (previously rotated) Greenwich meridian.


Except note that GRIB-2 uses southern pole of projection. Not sure how
that affects the algorithm.

The CF parameter north_pole_grid_longitude likely should be Angle of
rotation.



Anyone using this or have some insight on it ??




___
CF-metadata mailing list
CF-metadata@cgd.ucar.edu
http://mailman.cgd.ucar.edu/mailman/listinfo/cf-metadata





___
CF-metadata mailing list
CF-metadata@cgd.ucar.edu
http://mailman.cgd.ucar.edu/mailman/listinfo/cf-metadata


Re: [CF-metadata] Rotated Pole projection possibly wrong

2013-03-06 Thread David Hassell
Hello,

It is a bit of mess! As I understand it, the full rotation described
is a sequence of rotations about three different axes:

   Z  = 90S - 90N
   Y' = 90W' - 90E'
   X = 0E - 180E

where it is understood that the definitions of the N-S and W-E axes
change after each rotation (hence the primes and double
primes). Therefore the order in which they are done matters.

I suspect that the usual and assumed order is Z, Y', X?

From the GRIB-2 stuff John posted the north_pole_grid_longitude gives
the rotation about the X axis.

Something like John's Angle of rotation seems right to me. 

The apparent lack of consistency between the parameter names irks
me. Perhaps one solution could be to:

i) Add some text to the conventions state the order (Z, Y', X, say)
   and direction of rotations.

ii) add three new, consistent, self-describing parameters (e.g.)

angle_of_rotation_z
angle_of_rotation_y
angle_of_rotation_x

iii) allow aliases for backwards compatibility

grid_north_pole_longitude = angle_of_rotation_z
grid_north_pole_latitude  = angle_of_rotation_y
north_pole_grid_longitude = angle_of_rotation_x   

Or would that just obfuscate things even more?

All the best,

David

P.S. If you have a copy available, there are some nice descriptions in
 Coordinate Systems and Map Projections by D. H. Maling

--
David Hassell
National Centre for Atmospheric Science (NCAS)
Department of Meteorology, University of Reading,
Earley Gate, PO Box 243,
Reading RG6 6BB, U.K.

Tel   : 0118 3785613
Fax   : 0118 3788316
E-mail: d.c.hass...@reading.ac.uk
___
CF-metadata mailing list
CF-metadata@cgd.ucar.edu
http://mailman.cgd.ucar.edu/mailman/listinfo/cf-metadata


Re: [CF-metadata] Rotated Pole projection possibly wrong

2013-03-06 Thread John Caron

Hi Burkhardt:

Of course we cant change the existing conventions as there are files 
that already use it. But we need to document it so that we know when we 
are doing it right.


Im cleaning up my code and noticed that i have two implementations, one 
for grib and one for CF. So Id like to merge the two.


John

On 3/6/2013 12:32 PM, Burkhardt Rockel wrote:

John,

the rotated pole definition was requested by me in 2002 (see the thread 
projections/default cell_methods). In this thread the present description was 
chosen.
Please note: the rotated pole grid mapping is NOT a PROJECTION!  it is a 
TRANSFORMATION (therefore the GRIB description is misleading). This is 
probably the reason why you did not find it in Proj4.
The a, b, c in the GRIB definition description is fine. Whether  the rotated 
South Pole (grib) or North Pole (CF-Conventions) is used in the definition does 
not matter, the values can be easily converted to each other.

The rotated pole grid mapping as defined in the CF-Conventions document is 
already used widely as netCDF output standard for regional climate models in 
several international projects. It is also the standard in the IPCC-AR5 CORDEX 
model output for regional climate models. Therefore I am not in favor of 
changing the definition. However, a more precise description in the 
CF-Conventions would certainly be helpful for those not familiar with the 
rotated pole transformation.

Regards
Burkhardt

---
Dr. Burkhardt Rockel
Helmholtz-Zentrum Geesthacht
Institute of Coastal Research / Group Regional Atmospheric Modeling
Max-Planck-Strasse 1
D-21502 Geesthacht
Germany
Phone: +49 4152 87 1803
Fax: +49 4152 87 4 1803
Email: Burkhardt.Rockel (at) hzg.de
www: http://rockel.staff.coast.hzg.de
coordinates: 53.40579 N, 10.428647 E
---


Am 06.03.2013 um 17:09 schrieb John Caron ca...@unidata.ucar.edu:


Hi Heiko:

So

grid_north_pole_longitude = 
normalizeLongitude180(longitudeOfSouthernPoleInDegrees)
grid_north_pole_latitude = -1 * longitudeOfSouthernPoleInDegrees

?

where does one find documentation on proj4's ob_tran routine?

thanks,
John



On 3/6/2013 8:03 AM, Heiko Klein wrote:

Hi John,

all our models use rotated_latitute_longitude. And we have a consistent way of 
translating it between grib, proj and netcdf-java, i.e. publicly available at:
http://thredds.met.no/thredds/catalog/metno/proff4km/default/catalog.html


grid_mapping_name: rotated_latitude_longitude
grid_north_pole_longitude: 156.0
grid_north_pole_latitude: 23.5
earth_radius: 6367470.0

to proj:
proj4: +proj=ob_tran +o_proj=longlat +lon_0=-24 +o_lat_p=23.5 +a=6367470.0 
+no_defs

from grib:
lonRot = longitudeOfSouthernPoleInDegrees
latRot = latitudeOfSouthernPoleInDegrees

+proj=ob_tran +o_proj=longlat +lon_0=  normalizeLongitude180(lonRot)   
+o_lat_p=  (-1 * latRot);



This is the set of parameters which have been useful. I'm not sure if they are 
named in a useful way. proj4's ob_tran even allows for a o_lon_p parameter, but 
this doesn't seem to be used by grib or CF.

Best regards,

Heiko



On 2013-03-06 15:18, John Caron wrote:

The Rotated Pole projection here:

http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.6/cf-conventions.html#appendix-grid-mappings


Rotated pole

grid_mapping_name = rotated_latitude_longitude


Map parameters:

  *

grid_north_pole_latitude

  *

grid_north_pole_longitude

  *

north_pole_grid_longitude - This parameter is option (default is 0).

Map coordinates:

The rotated latitude and longitude coordinates are identified by
the standard_name attribute values
grid_latitude and grid_longitude respectively.



is problematic. It has no reference to Proj.4, and is not in Snyder.

Its likely the same as the GRIB-2 GDS Template 3.1:

Grid Definition Template 3.1: Rotated Latitude/longitude (or
equidistant cylindrical, or Plate Carrée)

Octet No. Contents

   15-72 Same as Grid Definition Template 3.0 (see Note 1)

   73-76 Latitude of the southern pole of projection

77-80 Longitude of the southern pole of
projection

   81-84 Angle of rotation of projection

   85-nnList of number of points along
each meridian or parallel (These octets are only present for
quasi-regular grids as described in Note 3)

Notes:


(2) Three parameters define a general latitude/longitude
coordinate system, formed by a general rotation of the sphere. One
choice for these parameters is:

(a) The geographic latitude in degrees of the southern pole of
the coordinate system, θp for example.

(b) The geographic longitude in degrees of the southern pole of
the coordinate system, λp for example.

(c) The angle of rotation in degrees about the new polar axis
(measured clockwise when looking from the southern to the northern pole)
of the coordinate system, assuming the new axis to have been obtained by
first rotating the sphere 

Re: [CF-metadata] Rotated Pole projection possibly wrong

2013-03-06 Thread John Caron

Hi Seth:

Thanks for the description, its really helpful.

Im just wondering how we document this in the CF convention, so that 
implementers have something to check against. Is there a reference 
implementation we can use ?


John

On 3/6/2013 1:17 PM, Seth McGinnis wrote:

So there are an infinite number of ways to decompose a 3-D orientation into
three orthogonal rotations.  Specifying the projection in terms of the final
location of the north pole and an optional rotation around that point means
that CF is actually agnostic about exactly how you get there.

That said, given that it's being specified in that way, the decomposition that
maps most straightforwardly to those parameters is a Z-Y-X rotation sequence
using an extrinsic (space-fixed) coordinate system:

First rotate by angle phi around the north pole in space-fixed coordinates.
Then rotate by angle theta around 0N 90E in space-fixed coordinates.
Finally, rotate by angle psi around 0N 0E in space-fixed coordinates.

The location of the grid's north pole in the original lat-lon coordinate system
is now at longitude = phi, colatitude = minus theta.

(I'm not sure how psi relates to north_pole_grid_longitude mathematically; it
defaults to zero and I've never seen it in use, so I think it's included mostly
for completeness. I agree that we may well want to change it to something more
obvious like an azimuth angle.)

Switching to the GRIB2 convention where you use the south pole instead of the
north should be just a matter of adding 180 degrees to phi and negating theta.

Cheers,

--Seth


On Wed, 6 Mar 2013 18:29:47 +
  David Hassell d.c.hass...@reading.ac.uk wrote:

Hello,

It is a bit of mess! As I understand it, the full rotation described
is a sequence of rotations about three different axes:

   Z  = 90S - 90N
   Y' = 90W' - 90E'
   X = 0E - 180E

where it is understood that the definitions of the N-S and W-E axes
change after each rotation (hence the primes and double
primes). Therefore the order in which they are done matters.

I suspect that the usual and assumed order is Z, Y', X?

From the GRIB-2 stuff John posted the north_pole_grid_longitude gives
the rotation about the X axis.

Something like John's Angle of rotation seems right to me.

The apparent lack of consistency between the parameter names irks
me. Perhaps one solution could be to:

i) Add some text to the conventions state the order (Z, Y', X, say)
   and direction of rotations.

ii) add three new, consistent, self-describing parameters (e.g.)

angle_of_rotation_z
angle_of_rotation_y
angle_of_rotation_x

iii) allow aliases for backwards compatibility

grid_north_pole_longitude = angle_of_rotation_z
grid_north_pole_latitude  = angle_of_rotation_y
north_pole_grid_longitude = angle_of_rotation_x

Or would that just obfuscate things even more?

All the best,

David

P.S. If you have a copy available, there are some nice descriptions in
 Coordinate Systems and Map Projections by D. H. Maling

--
David Hassell
National Centre for Atmospheric Science (NCAS)
Department of Meteorology, University of Reading,
Earley Gate, PO Box 243,
Reading RG6 6BB, U.K.

Tel   : 0118 3785613
Fax   : 0118 3788316
E-mail: d.c.hass...@reading.ac.uk
___
CF-metadata mailing list
CF-metadata@cgd.ucar.edu
http://mailman.cgd.ucar.edu/mailman/listinfo/cf-metadata

___
CF-metadata mailing list
CF-metadata@cgd.ucar.edu
http://mailman.cgd.ucar.edu/mailman/listinfo/cf-metadata


Re: [CF-metadata] Rotated Pole projection possibly wrong

2013-03-06 Thread Heiko Klein

Hi John,


grid_north_pole_latitude = -1 * longitudeOfSouthernPoleInDegrees
is ok



The normalize function is just used because grib is 0-360, while proj is 
-180 to 180. But you need to flip longitude when going from north to south:


grid_north_pole_longitude = norm180(longitudeOfSouthernPoleInDegrees + 180)


It's documented in the proj documentation 
ftp://ftp.remotesensing.org/proj/proj.4.3.I2.pdf Chapter 4 on page 19.


Heiko


On 2013-03-06 17:09, John Caron wrote:

Hi Heiko:

So

grid_north_pole_longitude =
normalizeLongitude180(longitudeOfSouthernPoleInDegrees)
grid_north_pole_latitude = -1 * longitudeOfSouthernPoleInDegrees

?

where does one find documentation on proj4's ob_tran routine?

thanks,
John



On 3/6/2013 8:03 AM, Heiko Klein wrote:

Hi John,

all our models use rotated_latitute_longitude. And we have a
consistent way of translating it between grib, proj and netcdf-java,
i.e. publicly available at:
http://thredds.met.no/thredds/catalog/metno/proff4km/default/catalog.html


grid_mapping_name: rotated_latitude_longitude
grid_north_pole_longitude: 156.0
grid_north_pole_latitude: 23.5
earth_radius: 6367470.0

to proj:
proj4: +proj=ob_tran +o_proj=longlat +lon_0=-24 +o_lat_p=23.5
+a=6367470.0 +no_defs

from grib:
lonRot = longitudeOfSouthernPoleInDegrees
latRot = latitudeOfSouthernPoleInDegrees

+proj=ob_tran +o_proj=longlat +lon_0= 
normalizeLongitude180(lonRot)   +o_lat_p=  (-1 * latRot);



This is the set of parameters which have been useful. I'm not sure if
they are named in a useful way. proj4's ob_tran even allows for a
o_lon_p parameter, but this doesn't seem to be used by grib or CF.

Best regards,

Heiko



On 2013-03-06 15:18, John Caron wrote:

The Rotated Pole projection here:

http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.6/cf-conventions.html#appendix-grid-mappings



Rotated pole

grid_mapping_name = rotated_latitude_longitude


Map parameters:

  *

grid_north_pole_latitude

  *

grid_north_pole_longitude

  *

north_pole_grid_longitude - This parameter is option (default
is 0).

Map coordinates:

The rotated latitude and longitude coordinates are identified by
the standard_name attribute values
grid_latitude and grid_longitude respectively.



is problematic. It has no reference to Proj.4, and is not in Snyder.

Its likely the same as the GRIB-2 GDS Template 3.1:

Grid Definition Template 3.1: Rotated Latitude/longitude (or
equidistant cylindrical, or Plate Carrée)

Octet No. Contents

   15-72 Same as Grid Definition Template 3.0 (see Note 1)

   73-76 Latitude of the southern pole of projection

77-80 Longitude of the southern pole of
projection

   81-84 Angle of rotation of projection

   85-nnList of number of points along
each meridian or parallel (These octets are only present for
quasi-regular grids as described in Note 3)

Notes:


(2) Three parameters define a general latitude/longitude
coordinate system, formed by a general rotation of the sphere. One
choice for these parameters is:

(a) The geographic latitude in degrees of the southern pole of
the coordinate system, θp for example.

(b) The geographic longitude in degrees of the southern pole of
the coordinate system, λp for example.

(c) The angle of rotation in degrees about the new polar axis
(measured clockwise when looking from the southern to the northern pole)
of the coordinate system, assuming the new axis to have been obtained by
first rotating the sphere through λp degrees about the geographic polar
axis, and then rotating through (90 + θp) degrees so that the southern
pole moved along the (previously rotated) Greenwich meridian.


Except note that GRIB-2 uses southern pole of projection. Not sure how
that affects the algorithm.

The CF parameter north_pole_grid_longitude likely should be Angle of
rotation.



Anyone using this or have some insight on it ??




___
CF-metadata mailing list
CF-metadata@cgd.ucar.edu
http://mailman.cgd.ucar.edu/mailman/listinfo/cf-metadata







--
Dr. Heiko Klein  Tel. + 47 22 96 32 58
Development Section / IT Department  Fax. + 47 22 69 63 55
Norwegian Meteorological Institute   http://www.met.no
P.O. Box 43 Blindern  0313 Oslo NORWAY

___
CF-metadata mailing list
CF-metadata@cgd.ucar.edu
http://mailman.cgd.ucar.edu/mailman/listinfo/cf-metadata