[R-sig-Geo] Raster projection alignment issues
Thomas Adams
tea3rd at gmail.com
Sat Dec 20 18:21:24 CET 2014
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 <
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.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
> 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
>
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list