Coordinate System – Convert Over 360 Degree Range to +/-180 Longitude in R

coordinate systemcoordinatesr

I'm working with data on sea surface currents from ftp://podaac-ftp.jpl.nasa.gov/allData/oscar/preview/L4/oscar_third_deg The data has the latitudes and longitudes saved in the following format:

latitude  Size:481
units: degrees-north
long_name: Latitude
longitude  Size:1201
units: degrees-east
long_name: Longitude

head(lon)
[1] 20.00000 20.33333 20.66667 21.00000 21.33333 21.66667
tail(lon)
[1] 418.3333 418.6667 419.0000 419.3333 419.6667 420.0000

How can I convert this into latitude-longitude units in R?

The lon object goes from 20 to 420 whereas I want it to be from -180 to 180.

I have seen Convert North and East coordinates to longitude and latitude in R at Stack Overflow but I am unable to adapt it to my setting.

x <- nc_open("oscar_vel1992.nc")
# lon
lon <- ncvar_get(x,"longitude")
#lat
lat <- ncvar_get(x,"latitude")
z <-as.numeric(char2dms(lat))

Best Answer

The data goes from 20 to 420, which is including some areas more than once because the metadata says:

NC_GLOBAL#NOTE2=Longitude extends from 20 E to 420 E to avoid a break in major ocean basins. Data repeats in overlap region.

So we need to split the data into the 20 to 360 section and the 360 to 380 section and put that 360 to 380 section on the left of the 20 to 360 section. We don't need 380 to 420 because that's overlap.

So...

library(raster)

define two extents to cut the raster with. The "Y" extent just has to be bigger than the raster so I go up to 90:

e1 = extent(c(xmin=20,xmax=360,ymin=-90,ymax=90))
e2 = extent(c(xmin=360, xmax=360+20, ymin=-90, ymax=90))

Load the raster - I'm loading just the "vm" layer, you'll have to repeat over the others:

d = raster("./oscar_vel7001.nc",varname="vm")

Now split the bits that we want using the extents:

d1 = crop(d,e1)
d2 = crop(d,e2)

Next shift the bit on the right to fill the space from 0 to 20:

xmin(d2)=xmin(d2)-360
xmax(d2)=xmax(d2)-360

merge will paste the two parts together, giving us a raster from 0 to 360, and rotate will shift it to be -180 to 180:

dm = rotate(merge(d1,d2))
plot(d)
plot(dm)

shows the original shifted overlapping map and the corrected -180 to +180 map:

Related Solutions

Latitude and Longitude Transformation – How to Manually Convert Rotated Lat/Lon to Regular Lat/Lon

Manually reversing the rotation should do the trick; there should be a formula for rotating spherical coordinate systems somewhere, but since I can't find it, here's the derivation ( ' marks the rotated coordinate system; normal geographic coordinates use plain symbols):

First convert the data in the second dataset from spherical (lon', lat') to (x',y',z') using:

x' = cos(lon')*cos(lat')
y' = sin(lon')*cos(lat')
z' = sin(lat')

Then use two rotation matrices to rotate the second coordinate system so that it coincides with the first 'normal' one. We'll be rotating the coordinate axes, so we can use the axis rotation matrices. We need to reverse the sign in the ϑ matrix to match the rotation sense used in the ECMWF definition, which seems to be different from the standard positive direction.

Since we're undoing the rotation described in the definition of the coordinate system, we first rotate by ϑ = -(90 + lat0) = -55 degrees around the y' axis (along the rotated Greenwich meridian) and then by φ = -lon0 = +15 degrees around the z axis):

x   ( cos(φ), sin(φ), 0) (  cos(ϑ), 0, sin(ϑ)) (x')
y = (-sin(φ), cos(φ), 0).(  0     , 1, 0     ).(y')
z   ( 0     , 0     , 1) ( -sin(ϑ), 0, cos(ϑ)) (z')

Expanded, this becomes:

x = cos(ϑ) cos(φ) x' + sin(φ) y' + sin(ϑ) cos(φ) z'
y = -cos(ϑ) sin(φ) x' + cos(φ) y' - sin(ϑ) sin(φ) z'
z = -sin(ϑ) x' + cos(ϑ) z'

Then convert back to 'normal' (lat,lon) using

lat = arcsin(z)
lon = atan2(y, x)

If you don't have atan2, you can implement it yourself by using atan(y/x) and examining the signs of x and y

Make sure that you convert all angles to radians before using the trigonometric functions, or you'll get weird results; convert back to degrees in the end if that's what you prefer...

Example (rotated sphere coordinates ==> standard geographic coordinates):

southern pole of the rotated CS is (lat0, lon0)

(-90°, *) ==> (-35°, -15°)
prime meridian of the rotated CS is the -15° meridian in geographic (rotated 55° towards north)

(0°, 0°) ==> (55°, -15°)
symmetry requires that both equators intersect at 90°/-90° in the new CS, or 75°/-105° in geographic coordinates

(0°, 90°) ==> (0°, 75°)
(0°, -90°) ==> (0°,-105°)

EDIT: Rewritten the answer thanks to very constructive comment by whuber: the matrices and the expansion are now in sync, using proper signs for the rotation parameters; added reference to the definition of the matrices; removed atan(y/x) from the answer; added examples of conversion.

EDIT 2: It is possible to derive expressions for the same result without explicit tranformation into cartesian space. The x, y, z in the result can be substituted with their corresponding expressions, and the same can be repeated for x', y' and z'. After applying some trigonometric identities, the following single-step expressions emerge:

lat = arcsin(cos(ϑ) sin(lat') - cos(lon') sin(ϑ) cos(lat'))
lon = atan2(sin(lon'), tan(lat') sin(ϑ) + cos(lon') cos(ϑ)) - φ

[GIS] Given a set of coordinates, how can I calculate the minimum bounds

Hope i understand the question correctly...

We can solve the problem for longitude and latitude separately, so i will take your example with the longitudes: -76, -135 and 164.

First i would order them:

-135, -76, 164

Then i would add the left most coordinate to the right again: -135 + 360 = 225

-135, -76, 164, 225

Now we can calculate the gaps between the coordinates:

-135 (59) -76 (240) 164 (61) 225
             .......

The biggest gap (240) must be the boundary of the minimum bounding box, the part that does not belong to the box. The dotted line is the biggest part of the circle we can spare out. In our example that means, the boundary box starts with 164, includes -135 and ends with -76.

Best Answer

Related Solutions

Latitude and Longitude Transformation – How to Manually Convert Rotated Lat/Lon to Regular Lat/Lon

[GIS] Given a set of coordinates, how can I calculate the minimum bounds

Related Question