Align Raster with boundary shapefile: selecting right proj4 R

Question

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)

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

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/

Answer 1

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)