Raster projection alignment issues

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

Raster projection alignment issues

Wall, Wade A ERDC-RDE-CERL-IL
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.html#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
[hidden email]




        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Raster projection alignment issues

Robert Hijmans
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
<[hidden email]> 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.html#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
> [hidden email]
>
>
>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Raster projection alignment issues

Thomas Adams-2
In reply to this post by Wall, Wade A ERDC-RDE-CERL-IL
Wade,

I was able to import the data into GRASS GIS (which uses gdal) without any
problems (ultimately); the issue you're having is that the ACTUAL map
projection is either not specified -OR- is specified as 'unprojected',
which is what I think is the case. In reality, the map projection is
latitude-longitude and the data is written that way within the netcdf file.
In GRASS GIS I simply imported the data (a 1-commandline step) directly
into a latitude-longitude GRASS GIS LOCATION, which I already had setup. I
tested by over-laying U.S. state boundaries, which align perfectly. That
the crs is defined as Albers Equal Area in the data set is CLEARLY an error.

Personally, I think you would be far better served to use something like
GRASS GIS for this kind of work; I have had a lot of success doing this.
One very nice feature is that GRASS GIS plays together with R wonderfully,
through the spgrass6 contributed package. Going this direction, allows one
to read/write GRASS raster and vector files into/from R and do any and all
the statistical manipulations you want. Additionally, I think for many
things, you'll find you don't need R, like simple data visualization and
most raster arithmetic.

Best,
Tom


On Fri, Dec 19, 2014 at 7:34 AM, Wall, Wade A ERDC-RDE-CERL-IL <
[hidden email]> 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.html#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
> [hidden email]
>
>
>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: [EXTERNAL] Re: Raster projection alignment issues

Wall, Wade A ERDC-RDE-CERL-IL
In reply to this post by Robert Hijmans
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:[hidden email]]
Sent: Saturday, December 20, 2014 12:24 AM
To: Wall, Wade A ERDC-RDE-CERL-IL
Cc: [hidden email]
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 <[hidden email]> 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
> [hidden email]
>
>
>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo