Canadian census map division in R

I'm very new to R and to mapping, and I would like to create a map of certain data. I have a set of data named "D.Montreal" which shows which canadian census division visitors to Montreal came from in 2010. I would like to create a map using this data to show how many people came from different regions, maybe by coloring the map in different colors depending on the number of people.

Is there a package I can use to have the delimitations of census divisions in Canada?

So far, I've made the outline of Canada :

map("worldHires","Canada", xlim=c(-141,-53), ylim=c(40,85), col="grey90", fill=TRUE)

Thank you so much!

Definitely a broad question and working through the spatial task view (as Roman suggested) is a great starting point. I prefer working with ggplot when it comes to mapping, though it takes some elbow grease. Some of the following won't really make sense until you hit some of the spatial task topics, but it essentially grabs a shapefile from Statistics Canada, simplifies it (so the polygons don't take forever to load/plot) and then colors them in by area:

library(rgeos)
library(rgdal)
library(maptools)
library(sp)
library(ggplot2)

# decent, uncluttered map theme (needs devtools package tho)
devtools::source_gist("https://gist.github.com/hrbrmstr/33baa3a79c5cfef0f6df")

# grab the file from "Statistics Canada"
download.file("http://www12.statcan.gc.ca/census-recensement/2011/geo/bound-limit/files-fichiers/gcd_000b11a_e.zip",
              destfile="gcd_000b11a_e.zip")

unzip("gcd_000b11a_e.zip")

# this simplifies the polygons so they load/plot faster
system("ogr2ogr canada.shp gcd_000b11a_e.shp -simplify 0.01")

# what layers do we have? you can use this to check
#   ogrListLayers("gcd_000b11a_e/canada.shp")
# but there are none, so the shapefile is the layer

canada <- readOGR("gcd_000b11a_e/","canada")

# do this to see what's available from an "identifier" standpoint
# "CDNAME" seems to be the census district name
# "PRNAME" seems to be the province name
# str(canada@data)

# rig up  some data
# make a data frame of census division areas
# you can assign as many value columns as you like
# they get merged in later and can be used as the fill level
# we'll use the area as the fill level
map_areas <- data.frame(id=canada@data$CDNAME,
                        area=sapply(slot(canada, "polygons"), slot, "area") )

# this takes a while, but it makes a data frame for use with
# ggplot and lets us use the census division name for doing things
# like applying colors
canada_map <- fortify(canada, region="CDNAME")

# merge in areas
canada_map <- merge(canada_map, map_areas, by="id")

gg <- ggplot()
gg <- gg + geom_map(data=canada_map, map=canada_map,
                    aes(map_id=id, x=long, y=lat, group=group, fill=log1p(area)),
                    color="white", size=0.1)
gg <- gg + coord_map() # can choose other projections
gg <- gg + theme_map()
gg

To work with something like population data is very similar. I found some data for province population and the following makes a choropleth for those (which may be over-simplified, but I was aiming for making the census districts easier to plot) based on population.

province_pop <- data.frame(
  id=c("Newfoundland and Labrador / Terre-Neuve-et-Labrador",
       "Prince Edward Island / Île-du-Prince-Édouard",
       "Nova Scotia / Nouvelle-Écosse", "New Brunswick / Nouveau-Brunswick",
       "Quebec / Québec", "Ontario", "Manitoba", "Saskatchewan",
       "Alberta", "British Columbia / Colombie-Britannique", "Yukon",
       "Northwest Territories / Territoires du Nord-Ouest", "Nunavut"),
  population=c(526977.0, 146283.0, 942668.0, 753914.0, 8214672.0, 13678740.0, 
               1282043.0, 1125410.0, 4121692.0, 4631302.0, 36510.0, 43623.0, 36585.0))

canada_map <- fortify(canada, region="PRNAME")
canada_map <- merge(canada_map, province_pop, by="id")

gg <- ggplot()
gg <- gg + geom_map(data=canada_map, map=canada_map,
                    aes(map_id=id, x=long, y=lat, group=group, fill=population),
                    color="white", size=0.5)
gg <- gg + scale_fill_continuous(low="#ccebc5", high="#084081")
gg <- gg + coord_map()
gg <- gg + theme_map()
gg

You'll definitely want to take a look at the values of the columns in canada@data since, as you can see with the province names, they may contain both french and english versions.

I can fake some census district population data to show it's the same methodology:

fake_census_pop <- data.frame(
  id=unique(as.character(canada@data$CDNAME)),
  fake_pop=sample(500000:30000000, length(unique(as.character(canada@data$CDNAME)))))

canada_map <- fortify(canada, region="CDNAME")
canada_map <- merge(canada_map, fake_census_pop, by="id")

gg <- ggplot()
gg <- gg + geom_map(data=canada_map, map=canada_map,
                    aes(map_id=id, x=long, y=lat, group=group, fill=fake_pop),
                    color="white", size=0.1)
gg <- gg + scale_fill_continuous(low="#ccebc5", high="#084081")
gg <- gg + coord_map()
gg <- gg + theme_map()
gg