[GIS] projection differences in R

coordinate systemr

I am trying to reproject some climate data that I downloaded 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 try to reproject rast1 using the other objects (rast2), the extent does not change and even though the reprojected rast1 is in the same crs, the extents don't line up.

rast_new = projectRaster(rast1, rast2)

So I had 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, even though rast_new is in the same projection as rast2

+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

the plotting isn't the same. See image plots of two images with same crs

Just wondering what I can do to correct this? Needless to say, I am relatively new to understanding projections in R ( and projections in general!)

I have exported the two rasters and viewed in ArcGIS, to see if it was a plot() issue in R, but the images are the same as above, yet the projections are the same (according to ArcGIS 10.2). Even when I reproject the "bad" raster in ArcGIS, the projection does not appear to be in Albers Equal Area. When I import directly from the netCDF file into ArcGIS, the projection looks correct.

Best Answer

It looks like the first image is using the WSG84 geodetic reference system, whereas the second image is using GRS80, which is why is looks curved. This suggests the objects haven't been converted to the same CRS.

Can you provide the climate data, or a sample, so I or another user can look into this further?

Related Solutions

[GIS] Extract values from raster using point data with buffer in R

If buffer is in meters the code works, but if I use degrees, it won't.

Indeed; the buffer cannot be specified in degrees. ?extract, buffer says: "If the data are not projected (latitude/longitude), the unit should be meters. Otherwise it should be in map-units (typically also meters)."

Do you really want the buffer to be in degrees?

Is so, temporarily set the CRS to UTM, do the extract, and set it back.

cr <- crs(r)
crs(r) <- '+proj=utm +zone=10 +datum=WGS84'
e <- extract(r, xy, buffer=0.5)
crs(r) <- cr

Process Vector to Raster Faster with R – Speed Optimization Techniques

I tried to "parallelize" the function rasterize using the R package parallel in this way:

split the SpatialPolygonsDataFrame object in n parts
rasterize every part separately
merge all the parts into one raster

In my computer, the parallelized rasterize function took 2.75 times less than the no-parallelized rasterize function.

Note: the code below download a polygon shapefile (~26.2 MB) from the web. You can use any SpatialPolygonDataFrame object. This is only an example.

Load libraries and example data:

# Load libraries
library('raster')
library('rgdal')

# Load a SpatialPolygonsDataFrame example
# Load Brazil administrative level 2 shapefile
BRA_adm2 <- raster::getData(country = "BRA", level = 2)

# Convert NAMES level 2 to factor 
BRA_adm2$NAME_2 <- as.factor(BRA_adm2$NAME_2)

# Plot BRA_adm2
plot(BRA_adm2)
box()

# Define RasterLayer object
r.raster <- raster()

# Define raster extent
extent(r.raster) <- extent(BRA_adm2)

# Define pixel size
res(r.raster) <- 0.1

Figure 1: Brazil SpatialPolygonsDataFrame plot

Simple thread example

# Simple thread -----------------------------------------------------------

# Rasterize
system.time(BRA_adm2.r <- rasterize(BRA_adm2, r.raster, 'NAME_2'))

Time in my laptop:

# Output:
# user  system elapsed 
# 23.883    0.010   23.891

Multithread thread example

# Multithread -------------------------------------------------------------

# Load 'parallel' package for support Parallel computation in R
library('parallel')

# Calculate the number of cores
no_cores <- detectCores() - 1

# Number of polygons features in SPDF
features <- 1:nrow(BRA_adm2[,])

# Split features in n parts
n <- 50
parts <- split(features, cut(features, n))

# Initiate cluster (after loading all the necessary object to R environment: BRA_adm2, parts, r.raster, n)
cl <- makeCluster(no_cores, type = "FORK")
print(cl)

# Parallelize rasterize function
system.time(rParts <- parLapply(cl = cl, X = 1:n, fun = function(x) rasterize(BRA_adm2[parts[[x]],], r.raster, 'NAME_2')))

# Finish
stopCluster(cl)

# Merge all raster parts
rMerge <- do.call(merge, rParts)

# Plot raster
plot(rMerge)

Figure 2: Brazil Raster plot