[R-sig-Geo] projectRaster across the antimeridian

[Bacou, Melanie]

Tony,
I can confirm this seems like an issue with `projectRaster()`. Using 
`gdalwarp` instead on the same tile seems to work.

# Download geoTIFF
BS.terra <- 
curl::curl_download("https://lance.modaps.eosdis.nasa.gov/imagery/subsets/?subset=BeringSea.2016182.terra.721.2km.tif",
   "BeringSea.2016182.terra.721.2km.tif")
BS.terra <- brick("BeringSea.2016182.terra.721.2km.tif")

# Reproject using `gdalwarp`
system("gdalwarp -overwrite -s_srs EPSG:4326 -t_srs EPSG:3035 -multi -of 
GTiff BeringSea.20e16182.terra.721.2km.tif 
BeringSea.2016182.terra.721.2km.azi.tif")

# Load reprojected raster
BS.terra.sa1 <- brick("BeringSea.2016182.terra.721.2km.azi.tif")

# Reproject with projectRaster()
BS.terra.sa2 <- projectRaster(from=BS.terra, crs=CRS("+init=epsg:3035"))

# Plot
par(mfrow=c(2,1))
plotRGB(BS.terra)
plotRGB(BS.terra.sa1)

# This one errors out
plotRGB(BS.terra.sa2)
 > Error in grDevices::rgb(RGB[, 1], RGB[, 2], RGB[, 3], alpha = alpha, 
max = scale) :
   color intensity -0.00130533, not in [0,1]

See https://dl.dropboxusercontent.com/u/30925475/Capture.PNG

--Mel.

On 7/1/2016 2:48 PM, Fischbach, Anthony wrote:
> The solution to the previous post addressing projection of rasters across
> the antimeridian
> (
> http://r-sig-geo.2731867.n2.nabble.com/coercing-marmap-bathy-object-to-raster-spanning-the-antimeridian-td7589828.html#a7589841
> )
> no-longer functions
> (as of 24 June 2016).
>
> Below I offer permutations of the solution offered in the previous post,
> to no avail.
>
> Is this a bug in the raster package, or may my code be adjusted to produce
> a projected raster that includes both the western and eastern hemisphere of
> the original image?
>
> -Tony
>
> require(sp)
> require(raster)
> setwd('D:/Temp')
> ## Define projections
> prj.geo<-CRS('+proj=longlat +datum=WGS84')
> prj.StudyArea <- CRS("+proj=aeqd +lat_0=70 +lon_0=-170")
> prj.StudyAreaPositive <- CRS("+proj=aeqd +lat_0=70 +lon_0=190")
> ## Download geoJPEG
> Yesterday <- Sys.Date() -1
> cat(URL<- format(Yesterday, '
> http://lance-modis.eosdis.nasa.gov/subsets/index.php?subset=BeringSea.%Y%j.terra.721.2km.jpg'),
> fill=T)
> cat(DestFile<-format(Yesterday,'BeringSea.%Y%j.terra.721.2km.jpg' ,
> sep=""), fill=T)##
> try(ds<-download.file(url=URL, destfile=DestFile, method="libcurl", quiet =
> F, mode = "wb", cacheOK = F))
> ## Cast as raster stack
> BS.terra<-stack(DestFile)
> ## Assign projection
> projection(BS.terra)<- prj.geo
> ymax(BS.terra) <- 70
> ymin(BS.terra) <- 58
> xmax(BS.terra) = 360-155
> xmin(BS.terra) = 360-190
> ## project the raster stack
> (BS.terra.sa1<-projectRaster(from=BS.terra,  res=5000,
> crs=prj.StudyAreaPositive, method='ngb', over=TRUE))
> (BS.terra.sa2<-projectRaster(from=BS.terra,  res=5000,
> crs=prj.StudyAreaPositive, method='ngb', over=F))
> (BS.terra.sa3<-projectRaster(from=BS.terra,  res=5000, crs=prj.StudyArea,
> method='ngb', over=TRUE))
> (BS.terra.sa4<-projectRaster(from=BS.terra,  res=5000, crs=prj.StudyArea,
> method='ngb', over=F))
> ## plot
> par(mfrow=c(5,1))
> plotRGB(BS.terra)
> plotRGB(BS.terra.sa1)
> plotRGB(BS.terra.sa2)
> plotRGB(BS.terra.sa3)
> plotRGB(BS.terra.sa4)
> #
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo