1
library(maps)
library(mapproj)
library(mapdata)
library(rgeos)
library(maptools)
library(sp)
library(raster)
library(rgdal)

Drawing inspiration from here https://stackoverflow.com/questions/10763421/r-creating-a-map-of-selected-canadian-provinces-and-u-s-states, I would like to align the raster and shapefile below. Information from the data providers say datis on a polar stereographicprojection. However, I could not get it to my shapefile to align with the rasterLayer. Here is my attempt so far:

dat file~1 KB

#saveRDS(dat1,'dat.rds')

dat<- readRDS('dat.rds')

# can0<-getData('GADM', country="CAN", level=0) # Canada
can1<-getData('GADM', country="CAN", level=1) # provinces

## Specify a geographic extent for the map
## by defining the top-left and bottom-right geographic coordinates
mapExtent <- rbind(c(-156, 80), c(-68, 19))

## Specify the required projection using a proj4 string
## Use http://www.spatialreference.org/ to find the required string
## Polyconic for North America
newProj <- CRS("+proj=poly +lat_0=0 +lon_0=-100 +x_0=0 
            +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")

## Project the map extent (first need to specify that it is longlat) 
mapExtentPr <- spTransform(SpatialPoints(mapExtent, 
                                         proj4string=CRS("+proj=longlat")),
                           newProj)

## Project other layers
can1Pr <- spTransform(can1, newProj)
#us1Pr <- spTransform(us1, newProj) 

## Plot each projected layer, beginning with the projected extent
extent(dat)=extent(can1Pr)
plot(mapExtentPr, pch=NA)
plot(dat, add=TRUE)
plot(can1Pr, border="black", col=NA, add=TRUE)

#===================

enter image description here

Additional information about the data:

The CANGRD grid is in polar stereographic projection with a 50 km spatial resolution. The grid is a 125 (columns) by 95 (rows) matrix, where the SW corner (0,0) is at 40.0451°N latitude and 129.8530°W longitude. The projection is true at 60.0°N and centered on 110.0°W. The file ‘CANGRD_points_LL.txt’ lists the latitudes and longitudes for each grid point.

The general format of the ‘YYYYDD.grd’ file is:
Id – ‘DSAA’ identifies the file as an ASCII grid file
nx ny - nx is the integer number of grid lines along the X axis (columns)
        ny is the integer number of grid lines along the Y axis (rows)
xlo xhi - xlo is the minimum X value of the grid
          xhi is the maximum X value of the grid
ylo yhi - ylo is the minimum Y value of the grid
          yhi is the maximum Y value of the grid
zlo zhi - are the minimum and maximum Z values of the grid. 

data source: ftp://ccrp.tor.ec.gc.ca/pub/EC_data/CANGRD/

2

EDIT --- Here's the solution using the additional information you provided

library(raster)    
library(proj4)

# Define projection based on CANGRD description
proj4s <- "+proj=stere +lat_0=90 +lat_ts=60 +lon_0=-110 +x_0=0 +y_0=0 +datum=WGS84"
sw_ll <- c(-129.853, 40.0451)

# convert sw corner to projected coordinates
sw_proj <- project(sw_ll, proj4s)

# Define extent
xmin <- sw_proj[1]
ymin <- sw_proj[2]
xmax <- xmin + 125 * 50000
ymax <- ymin + 95 * 50000
e <- extent(xmin, xmax, ymin, ymax)

# Read raster and assign proper CRS and extent
dat <- raster('/path/to/dat.nc')
crs(dat) <- proj4s
extent(dat) <- e

# Get province borders and project it to same CRS than raster
can1 <- getData('GADM', country="CAN", level=1)
can_proj <- spTransform(can1, CRS(projection(dat)))

# Plot result
plot(dat)
plot(can_proj, border="black", col=NA, add=TRUE)

Original answer

R thinks that your raster file is in longlat. If you're certain that it's in polar stereographic, the following should work.

dat <- readRDS('dat.rds')

crs(dat) <- "+proj=poly +lat_0=0 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"

can_ll <- getData('GADM', country="CAN", level=1)

can_proj <- spTransform(can_ll, CRS(projection(dat)))

plot(dat)
plot(can_proj, border="black", col=NA, add=TRUE)
  • did you get the shapefile to align with raster? Your solution does not yet get it to align. For example, extent (dat) : -141, -52, 41, 83.5 (xmin, xmax, ymin, ymax) and extent (can_proj) : -2123593, 3383194, 4760230, 9356298 (xmin, xmax, ymin, ymax) are different so they cannot align. Setting extent(dat)=extent(can_proj) doesn't work either. Something might be wrong in my code? – code123 Apr 26 '17 at 22:01
  • I couldn't test it since the rds file you provided only contains the metadata of the raster file; run readAll() before saveRDS() if you want to save data too. – Loïc Dutrieux Apr 26 '17 at 22:22
  • The extent values for dat are definitely in degrees. Try passing crs="+proj=poly +lat_0=0 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" to raster() when reading the file from the netcdf archive. – Loïc Dutrieux Apr 26 '17 at 22:28
  • here is a link to dat data in .nc format. dropbox.com/s/a6bm2selasly584/dat.nc?dl=0 – code123 Apr 26 '17 at 22:30
  • The problem is with that netcdf file that has a longlat CRS but whose grid is clearly not longlat. Your best chance is probably to manually adjust the extent of your raster. If you know the top left corner coordinate of your raster and its resolution that should be doable. – Loïc Dutrieux Apr 26 '17 at 22:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.