make rotate generic, add sf & sfc methods

DESCRIPTION

@@ -109,6 +109,7 @@ Collate:
     'proxy.R'
     'factors.R'
     'rasterize.R'
+    'rotate.R'
     'subset.R'
     'warp.R'
     'aggregate.R'

NAMESPACE

@@ -132,6 +132,9 @@ S3method(st_normalize,stars_proxy)
 S3method(st_redimension,nc_proxy)
 S3method(st_redimension,stars)
 S3method(st_redimension,stars_proxy)
+S3method(st_rotate,sf)
+S3method(st_rotate,sfc)
+S3method(st_rotate,stars)
 S3method(st_sample,stars)
 S3method(st_sample,stars_proxy)
 S3method(st_set_bbox,dimensions)

R/rotate.R

@@ -0,0 +1,96 @@
+toRad = function(x) x * pi / 180
+
+toDeg = function(x) x * 180 / pi
+
+rotate = function(lonlat, lon0, lat0, north = TRUE) {
+	# https://gis.stackexchange.com/questions/10808/manually-transforming-rotated-lat-lon-to-regular-lat-lon/14445 , by Miha!
+	lon0 = toRad(lon0)
+	lat0 = toRad(lat0)
+	if (north) {
+		lat0 = -lat0
+		lon0 = pi + lon0
+	}
+	lonlat = toRad(lonlat)
+	lon = lonlat[,1]
+	lat = lonlat[,2]
+	vartheta = -(pi/2 + lat0)
+	varphi = -lon0
+	toDeg(cbind(
+		atan2(sin(lon), tan(lat) * sin(vartheta) + 
+				cos(lon) * cos(vartheta)) - varphi,
+		asin(cos(vartheta) * sin(lat) - cos(lon) * sin(vartheta) * cos(lat))
+	))
+}
+
+#' @export
+st_rotate = function(.x, lon0, lat0, north, ...) UseMethod("st_rotate")
+
+
+#' Transform rotated pole long/lat regular grid to unrotated curvilinear grid
+#'
+#' Transform rotated long/lat regular grid to unrotated curvilinear grid
+#' @param .x object of class \code{stars}
+#' @param lon0 longitude of the rotated pole in degrees
+#' @param lat0 latitude of the rotated pole in degrees
+#' @param north logical; if \code{TRUE} the pole refers to the North pole, otherwise the South pole
+#' @param ... ignored
+#' @returns curvilinear stars object with coordinates in regular long/lat (North pole at lat=90)
+#' @export
+#' @name st_rotate
+#' @examples
+#' if (require("starsdata") && require("maps")) {
+#'   # data downloaded from https://esgf-data.dkrz.de/search/cosmo-rea/
+#'   nc = "netcdf/ts_EUR-6km_ECMWF-ERAINT_REA6_r1i1p1f1_COSMO_v1_mon_201801-201812.nc"
+#'   f = system.file(nc, package = "starsdata")
+#'   m = read_mdim(f, "ts")
+#'   print(m)
+#'   # NOTE this function is obsolete when reading m as
+#'   # m = read_mdim(f, "ts", curvilinear = c("longitude", "latitude"))
+#'   if (require(RNetCDF)) {
+#'      x = open.nc(f)
+#'      lon = att.get.nc(x, "rotated_latitude_longitude", "grid_north_pole_longitude")
+#'      lat = att.get.nc(x, "rotated_latitude_longitude", "grid_north_pole_latitude")
+#'      close.nc(x)
+#'      print(c(lon = lon, lat = lat))
+#'   } else {
+#'      lon = -162
+#'      lat = 39.25
+#'   } 
+#'   m1 = st_rotate(m, lon, lat)
+#'   print(m1)
+#'   h = function() maps::map(add = TRUE)
+#'   plot(m1, downsample = c(10, 10, 5), axes = TRUE, hook = h, mfrow = c(1, 2)) 
+#'     # curvilinear grid: downsample for plotting speed
+#'   m2 = st_warp(m1, crs = st_crs("OGC:CRS84"), threshold = .1)
+#'   plot(m2, hook = h, mfrow = c(3, 4)) # regular grid: plots fast
+#' }
+st_rotate.stars = function(.x, lon0, lat0, north = TRUE, ...) {
+	stopifnot(is.na(st_crs(.x)) || st_is_longlat(.x), 
+			  !is_curvilinear(.x), 
+			  is.logical(north), 
+			  !is.na(north),
+			  lat0 <= 90,
+			  lat0 >= -90,
+			  length(list(...)) == 0)
+	if (has_sfc(.x))
+		st_set_dimensions(.x, which_sfc(.x), values = rotate(st_geometry(.x), lon0, lat0, north))
+	else {
+		.x = st_upfront(.x)
+		d = dim(.x)
+		ccr = rotate(st_coordinates(st_dimensions(.x)[1:2]), lon0, lat0, north)
+		st_as_stars(.x, curvilinear = setNames(list(matrix(ccr[,1], d[1], d[2]),
+					  matrix(ccr[,2], d[1], d[2])), names(d)[1:2]))
+	}
+}
+
+#' @export
+st_rotate.sfc = function(.x, lon0, lat0, north = TRUE, ...) {
+	r = rapply(.x, rotate, how = "replace", lon0 = lon0, lat0 = lat0, north = north)
+	st_set_crs(r, NA_crs_)
+}
+
+#' @export
+st_rotate.sf = function(.x, lon0, lat0, north = TRUE, ...) {
+	st_set_geometry(.x, st_rotate(st_geometry(.x), lon0, lat0, north))
+}
+

R/warp.R

@@ -252,71 +252,3 @@ st_warp = function(src, dest, ..., crs = NA_crs_, cellsize = NA_real_, segments
 	}
 	ret
 }
-
-rotate = function(lon, lat, lon0, lat0, north = TRUE) {
-	# https://gis.stackexchange.com/questions/10808/manually-transforming-rotated-lat-lon-to-regular-lat-lon/14445 , by Miha!
-	if (north) {
-		lat0 = -lat0
-		lon0 = pi + lon0
-	}
-	vartheta = -(pi/2 + lat0)
-	varphi = -lon0
-	cbind(
-		lon = atan2(sin(lon), tan(lat) * sin(vartheta) + 
-				cos(lon) * cos(vartheta)) - varphi,
-		lat = asin(cos(vartheta) * sin(lat) - cos(lon) * sin(vartheta) * cos(lat))
-	)
-}
-
-#' Transform rotated pole long/lat regular grid to unrotated curvilinear grid
-#'
-#' Transform rotated long/lat regular grid to unrotated curvilinear grid
-#' @param .x object of class \code{stars}
-#' @param lon0 longitude of the rotated pole in degrees
-#' @param lat0 latitude of the rotated pole in degrees
-#' @param north logical; if \code{TRUE} the pole refers to the North pole, otherwise the South pole
-#' @returns curvilinear stars object with coordinates in regular long/lat (North pole at lat=90)
-#' @export
-#' @examples
-#' if (require("starsdata") && require("maps")) {
-#'   nc = "netcdf/ts_EUR-6km_ECMWF-ERAINT_REA6_r1i1p1f1_COSMO_v1_mon_201801-201812.nc"
-#'   f = system.file(nc, package = "starsdata")
-#'   m = read_mdim(f, "ts")
-#'   print(m)
-#'   # NOTE this function is obsolete when reading m as
-#'   # m = read_mdim(f, "ts", curvilinear = c("longitude", "latitude"))
-#'   if (require(RNetCDF)) {
-#'      x = open.nc(f)
-#'      lon = att.get.nc(x, "rotated_latitude_longitude", "grid_north_pole_longitude")
-#'      lat = att.get.nc(x, "rotated_latitude_longitude", "grid_north_pole_latitude")
-#'      close.nc(x)
-#'      print(c(lon = lon, lat = lat))
-#'   } else {
-#'      lon = -162
-#'      lat = 39.25
-#'   } 
-#'   m1 = st_rotate(m, lon, lat)
-#'   print(m1)
-#'   h = function() maps::map(add = TRUE)
-#'   plot(m1, downsample = c(10, 10, 5), axes = TRUE, hook = h, mfrow = c(1, 2)) 
-#'     # curvilinear grid: downsample for plotting speed
-#'   m2 = st_warp(m1, crs = st_crs("OGC:CRS84"), threshold = .1)
-#'   plot(m2, hook = h, mfrow = c(3, 4)) # regular grid: plots fast
-#' }
-st_rotate = function(.x, lon0, lat0, north = TRUE) {
-	stopifnot(inherits(.x, "stars"), 
-			  is.na(st_crs(.x)) || st_is_longlat(.x), 
-			  !is_curvilinear(.x), 
-			  is.logical(north), 
-			  !is.na(north),
-			  lat0 <= 90,
-			  lat0 >= -90)
-	torad = function(x) x * pi / 180
-	todeg = function(x) x * 180 / pi
-	d = dim(.x)
-	n = prod(d[1:2])
-	cc = torad(st_coordinates(.x)[1:n, 1:2])
-	cc = todeg(rotate(cc[,1], cc[,2], torad(lon0), torad(lat0), north))
-	st_as_stars(.x, curvilinear = setNames(list(matrix(cc[,1], d[1], d[2]),
-					  matrix(cc[,2], d[1], d[2])), names(d)[1:2]))
-}