0

I would like to display the results in a raster times series object (rts) using levelplot function in rasterVis.

Here is a short code frome rts package:

    library(raster)
 library(rasterVis)
    library(rts)
    path <- system.file("external", package="rts")
    lst <- list.files(path=path,pattern='.asc$',full.names=TRUE)
    r <- stack(lst)
    d <- c("2000-02-01","2000-03-01","2000-04-01","2000-05-01") # corresponding dates to 4 rasters
    d <- as.Date(d)

    # creating a RasterStackTS object:
    rt <- rts(r,d)

To coerce from .rts to raster I used:

r=rt@raster



 proj=CRS("+proj=longlat +datum=NAD27 +no_defs +ellps=clrk66 +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat")

    proj4string(r) <- proj
wgs84.p4s <- CRS("+proj=longlat +datum=NAD83 +ellps=GRS80 +no_defs")

reprojecting r from UTM coordinates to latlong coordinates

rx <- projectRaster(from=r, crs=wgs84.p4s@projargs, method="ngb")

writeRaster(rt, file="myfile100.tif", format="GTiff", overwrite=TRUE)
ras = raster("myfile100.tif")

How can I read all 4 layers in myfile100.tif? ras = raster("myfile100.tif") reads only layer 1

proj4string(ras) <- proj
plot(ras)
levelplot(ras)

I get this error:

Error in projectExtent(from, projto) : cannot do this transformation
In addition: Warning message:
In rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[,  :
  146 projected point(s) not finite

How can I get around this? Maybe there is even a shorter and reasonable way to do this. Thanks for your suggestions.

www
  • 38,575
  • 12
  • 48
  • 84
code123
  • 2,082
  • 4
  • 30
  • 53
  • 1
    When reading in use `stack()` and it will read the tif in as a stack. Pondering the rest of your question. – Badger Jan 28 '16 at 17:21
  • @Badger ras = stack("myfile100.tif") reads all layers. Great! – code123 Jan 28 '16 at 17:46
  • Try `rx <- projectRaster(from=r, crs=wgs84.p4s, method="ngb")`, I've never used the @ part that you declared but its not to say that it is incorrect, I'm just curious of the behaviour. – Badger Jan 28 '16 at 17:52
  • @Badger I got the same error using `rx <- projectRaster(from=r, crs=wgs84.p4s, method="ngb")` – code123 Jan 28 '16 at 18:03
  • You should take a look at the `setZ` and `getZ` functions provided by the `raster` package. Maybe you could use them instead of the `rts` package. On the other hand, I think the error is due to the projection definition you are using. Are you sure that is the correct proj4 string? [This page](http://spatialreference.org/) will be useful to check it. – Oscar Perpiñán Jan 28 '16 at 23:32
  • @OscarPerpiñán I find `setZ` ,`zApply` and `getZ` very useful. But looks like these functions are still under development? How can I use `zApply` to sum `rasterStack` object with 30mins time step to daily time step? Then I can easily generate my maps using levelplot – code123 Jan 30 '16 at 02:05
  • You can use them safely. In my opinion, they are stable functions. I should suggest R. Hijmans to remove the warning about being under development. Regarding your question about the usage of `zApply`, please open a new question because comments are not the right place. Meanwhile, you have an example [here](https://stat.ethz.ch/pipermail/r-sig-geo/2013-April/018136.html). – Oscar Perpiñán Jan 30 '16 at 13:13

0 Answers0