Search code examples
rcoordinateszooutmproj

R: Lat-Lon to UTM for zoo time series


Given a zoo time series with columns "lat", "lon" and "value":

z <- zoo(...)
# anyone knows how to create a simple example series here?

I need to convert the lat and lon degrees to UTM in order to perform further operations like distance calculation. After a bit research I came up with a working snippet that converts everything to a data.frame:

d <- as.data.frame(z)
coordinates(d) <- ~ lon+lat
proj4string(d) <- CRS("+proj=longlat +datum=WGS84")
d.final <- spTransform(d, CRS(paste("+proj=utm +zone=",zone," ellps=WGS84",sep="")))

This removes the initial lat and lon columns and adds two new UTM-like lat and lon columns. Is there any simple way to do this directly with the zoo object?

Furthermore, this snippet clearly lacks information from the variable zone. From Determining UTM zone (to convert) from longitude/latitude we learn how to get the zone number from the longitude:

zone.number <- (floor((z$lon + 180)/6) %% 60) + 1

But I still do not now how to calculate the zone letter, eg. "N". How to obtain that? In general, I can expect the zone for all my time series to be constant.

This solution modified from an answer was what was used:

# coming from a zoo time series z with columns z$lat and z$lon
library(rgdal)
# assume that z stays within a zone
zone <- (floor((z$lon[1] + 180)/6) %% 60) + 1
# convert
utm <- project(merge(z$lon, z$lat), paste0("+proj=utm +zone=", zone))
# assign UTM values to new columns in z
z$utmx <- utm[,1]
z$utmy <- utm[,2]

Solution

  • Try this:

    library(zoo)
    library(rgdal)
    z <- zoo(cbind(lat = 12:14, lon = 10:12, value = 1:3)) # test data
    
    zone <- (floor((z$lon[1] + 180)/6) %% 60) + 1
    cbind(utm = project(z[, 2:1], paste0("+proj=utm +zone=", zone)), value = z[, 3])
    

    giving:

         utm.1   utm.2 value
    1 608864.2 1326751     1
    2 716922.6 1437988     2
    3 824104.0 1549780     3
    

    where utm.1 and utm.2 are x and y (or easting and northing) respectively. Note that the assumptions as to positive/negative are as in

    http://home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html

    and that the above R code output double checks to that page. It also double checks to the output in the link in the comments if you look at the standard UTM output on that page.

    Update: corrections and modified example.