[R-sig-Geo] [EXTERNAL] Re: Raster projection alignment issues

[Wall, Wade A ERDC-RDE-CERL-IL]

Thanks Dr. Hijmans,

That is a very simple solution. I had tried changing the extent by using other objects, but not directly manipulating the extent.

Wade 



-----Original Message-----
From: Robert J. Hijmans [mailto:r.hijmans at gmail.com] 
Sent: Saturday, December 20, 2014 12:24 AM
To: Wall, Wade A ERDC-RDE-CERL-IL
Cc: r-sig-geo at r-project.org
Subject: [EXTERNAL] Re: [R-sig-Geo] Raster projection alignment issues

Wade,

I do not think you are doing anything wrong. projectRaster gets confused because of the uncommon longitude coordinates used by climate models (0 to 360, rather than -180 to 180) You can use the 'rotate' function to fix that (ignore the warning), or by setting the extent, and then proceed. See below.

Robert

library(raster)
r <- raster(xmn=235.25, xmx=293, ymn=25.125, ymx=52.875, res=1/8, crs='+proj=longlat +datum=WGS84') rr <- rotate(r)

## or something like
## extent(r) <- extent(235.25-360, 293-360, 25.125, 52.875)

aea <- '+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0
+y_0=0 +datum=NAD83'

library(maptools)
data(wrld_simpl)
us <- wrld_simpl[wrld_simpl$ISO3=='USA', ]

library(rgdal)
usa <- spTransform(us, CRS(aea))
rra <- projectRaster(rr, crs=aea)

plot(rra)
plot(usa, add=TRUE)




by the way. this is almost certainly a bad idea:
values(rast_new) = values(rast1) ### This is necessary because all the values were NA

On Fri, Dec 19, 2014 at 6:34 AM, Wall, Wade A ERDC-RDE-CERL-IL <Wade.A.Wall at usace.army.mil> wrote:
> Hi all,
>
> I posted this over at gis.stackexchange.com, but thought this user group may be interested.
>
> I downloaded  raster data from  
> http://gdo-dcp.ucllnl.org/downscaled_cmip_projections/dcpInterface.htm
> l#Welcomein
>
> The BCSD-CMIP5 data is found under the tab "Projections: Complete Archives".
> The files are in netCDF format, but am running into some issues. The 
> original raster data (let's call it rast1) when imported has the 
> extent
>
> xmin 235.25
> xmax 293
> ymin 25.125
> ymax 52.875
>
> and a crs of +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 
> +towgs84=0,0,0
>
> The other objects (both sp and raster) I am working with contain 
> negative values for xmin and xmax
>
>
> xmin        : -122.8021
>
> xmax        : -75.38249
>
> ymin        : 31.78881
>
> ymax        : 47.15952
>
>
>
> When I attempt to reproject the raster data into Albers Equal Area 
> (+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 
> +datum=NAD83
>
> +units=m +no_defs +ellps=GRS80+towgs84=0,0,0)using the following code:
>
>
>
> rast_new = projectRaster(rast1, rast2) ## rast2 is a raster with the 
> crs for Albers Equal Area
>
>
>
> the extents still don't line up, so I have to do the following.
>
>
>
> rast_new = projectExtent(rast1,rast2)
>
> rast_new = projectRaster(rast1,rast_new)
>
> values(rast_new) = values(rast1) ### This is necessary because all the 
> values were NA
>
>
>
> However, rast_new, even though its crs is defined as Albers Equal Area, does not project correctly. Even when I export these files from R and open in ArcGIS, the issues persist. However, if I import the netCDF files directly into ArcGIS and reproject, there are no problems.
>
>
>
> Any clues as to what I am doing wrong?
>
>
>
> Thanks for any assistance.
>
>
> Wade A. Wall
> US Army ERDC-CERL
> P.O. Box 9005
> Champaign, IL  61826-9005
> 1-217-373-4420
> Wade.A.Wall at usace.army.mil
>
>
>
>
>         [[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