Showing posts with label Choropleth map. Show all posts
Showing posts with label Choropleth map. Show all posts

Oct 11, 2017

Random graphs (114): Line plot and choropleth

use ESS1-7e01, clear

// Country name variable
kountry cntry, from(iso2c) marker
rename NAMES_STD country

// Create No contract variable
gen nocontract = (wrkctra == 3) if !missing(wrkctra)

// Calculate per year and country
preserve
statsby cont = _b[_cons] contse = _se[_cons], by(country essround) clear: regress nocontract

replace cont = cont * 100             // Convert proportion into precentage
generate  lb = cont - (contse * 100)
generate  ub = cont + (contse * 100)

label define essround 1 "2002" 2 "2004" 3 "2006" 4 "2008" 5 "2010" 6 "2012" 7 "2014", modify

sort essround country
twoway (rarea lb ub essround) ///
       (connected cont essround), ///+
        by(country, legend(off) note(" " "{it:Source:} ESS 2002-2014", span) ///
        title("{bf:A}", justification(left) bexpand span)) ///
        xtitle("") xlabel(1/7, val ang(45)) ytitle("% no contract") ///
        xsize(6) ysize(6) ///
        name(figurea, replace) 
restore
  
// Calculate average % no contracts per country
collapse nocontract, by(cntry)
replace nocontract = nocontract * 100
format nocontract %6.1f
saveold cont, replace
  
// Read in map data
  // Source: shapefile from http://www.naturalearthdata.com/downloads/10m-cultural-vectors/
shp2dta using "maps\ne_10m_admin_0_countries", database(database)  coordinates(coordinates) genid(id) replace

  // Restrict coordinates to Europe
  // Source: https://en.wikipedia.org/wiki/Extreme_points_of_Europe
use coordinates, clear
replace _Y = . if _Y < 36
replace _Y = . if _Y > 71 & !missing(_Y)
replace _X = . if _X < -28 
replace _X = . if _X > 33 & !missing(_X)
saveold europecoordinates, replace

  // Create data set with Europe internal borders
clonevar id = _ID
merge m:1 id using database, nogenerate
keep if CONTINENT == "Europe"
saveold europe, replace

  // Merge map data with set with data to be plotted
use ISO_A2 NAME id using database, clear
rename ISO_A2 cntry
replace cntry = "FR" if NAME == "France" // Somehow not correct in map
replace cntry = "NO" if NAME == "Norway" // Somehow not correct in map
drop if cntry == "-99"
merge 1:1 cntry using cont, keep(matched) nogenerate

spmap nocontract using europecoordinates, id(id) xsize(6) ysize(6) ///
                                     polygon(data(europe)) legstyle(2) clnumber(9) ///
                                     legend(position (9) ring(0)) fcolor(Oranges) ///
                                     legorder(hilo) ocolor(none ..) ///
                                     title("{bf:B}", justification(left) bexpand span) ///
                                     legtitle("{bf:% no contract, 2002-2014 averages}") legjunction({c 150}) ///
                                     note(" " "{it:Source:} ESS 2002-2014", span) name(figureb, replace)
          
graph combine figurea figureb, row(1) xsize(12)

Jun 16, 2017

Random graphs (97): Choropleth map

// Calculate country data
use v9 weight_g c_abrv using ZA4800_v3-0-0.dta, clear
gen cntry = substr(c_abrv, 1, 2)
gen poorhealth = (inlist(v9, 4, 5) *100) if !missing(v9)
collapse poorhealth [pweight = weight_g], by(cntry)
label var poorhealth "Poor health"
format poorhealth %6.1f // Determine no. of decimal points
saveold poorhealth, replace

// Read in map data
  // Source: shapefile from http://www.naturalearthdata.com/downloads/10m-cultural-vectors/
shp2dta using "ne_10m_admin_0_countries", database(database)  coordinates(coordinates) genid(id) replace

  // Restrict coordinates to Europe
  // Source: https://en.wikipedia.org/wiki/Extreme_points_of_Europe
use coordinates, clear
replace _Y = . if _Y < 36
replace _Y = . if _Y > 71 & !missing(_Y)
replace _X = . if _X < -28 
replace _X = . if _X > 33 & !missing(_X)
saveold europecoordinates, replace

  // Create data set with Europe internal borders
clonevar id = _ID
merge m:1 id using database, nogenerate
keep if CONTINENT == "Europe"
saveold europe, replace

  // Merge map data with set with data to be plotted
use ISO_A2 NAME id using database, clear
rename ISO_A2 cntry
replace cntry = "FR" if NAME == "France" // Somehow not correct in map
replace cntry = "NO" if NAME == "Norway" // Somehow not correct in map
drop if cntry == "-99"
merge 1:1 cntry using poorhealth, keep(matched) nogenerate

spmap poorhealth using europecoordinates, id(id) xsize(6) ysize(4) ///
                                     polygon(data(europe)) legstyle(2) clnumber(5) fcolor(Terrain) ///
                                     legend(position (9) ring(0)) ///
                                     legorder(hilo) ocolor(none ..) ///
                                     legtitle("% poor health") legjunction({c 150})

Aug 20, 2012

Random graphs (2); Plotting country choropleths in R

 

 Steps in Stata

1) Create aggregate data set, for instance like this:
table cntry, c(mean sclmeet mean disc mean sclact) format (%9.4f)

2) Copy and paste resulting table to text editor or Excel to create a .csv file that can be used in R. Make sure that it's really comma-separated values, as in:
AT,4.6401,0.8263,2.7099
BE,4.8975,0.8292,2.5453
...
3) Also, add a header line in the first line of the file, so that it looks like:
Country,sclmeet,disc,sclact
AT,4.6401,0.8263,2.7099
BE,4.8975,0.8292,2.5453
...

4) Give it a name, for instance 'socialcontacts.csv,' and save it somewhere.

 

Steps in R

5) Open R. Install -rworldmaps-.
6) Set your working directory to where the data is and where you will save everything.
setwd("D:/YourDirectory/WhereTheDataIs")
7) Read in data and give it a name (here: "Meetsoc')
Meetsoc <-read.table("socialcontacts.csv", header=TRUE, sep=",")
attach(Meetsoc) # attach the data
names(Meetsoc) # give it names
Meetsoc[1:5,] # look at first 5 lines to see if it worked
8) Load 'rworldmap.'
library(rworldmap)
9) Join data to map and give it some name (here: 'MA1B')
MA1B <- joinCountryData2Map(Meetsoc, joinCode = "ISO2",
nameJoinColumn = "Country")
10) Produce map. xlim and ylim are the longitude and latitude. Missing country and ocean colors can be changed however one likes them.
# First map
par(mai = c(0, 0, 0.2, 0), xaxs = "i", yaxs = "i")

mapCountryData(MA1B, xlim=c(-15,29), ylim=c(40,70),
nameColumnToPlot = "sclmeet", catMethod = "categorical",
mapTitle = "Frequency of social meetings", colourPalette = "heat",
oceanCol = "lightblue", missingCountryCol = "white", borderCol =
"black", addLegend = F)

# Second map
par(mai = c(0, 0, 0.2, 0), xaxs = "i", yaxs = "i")

mapCountryData(MA1B, xlim=c(-15,29), ylim=c(40,70),
nameColumnToPlot = "disc", catMethod = "categorical",
mapTitle = "Prevalence of discussion partners",
colourPalette = "heat", oceanCol = "lightblue",
missingCountryCol = "white", borderCol = "black", addLegend = F)

# Third map
par(mai = c(0, 0, 0.2, 0), xaxs = "i", yaxs = "i")

mapCountryData(MA1B, xlim=c(-15,29), ylim=c(40,70),
nameColumnToPlot = "sclact", catMethod = "categorical",
mapTitle = "Frequency of social activities", 
colourPalette = "heat", oceanCol = "lightblue",
missingCountryCol = "white", borderCol = "black", addLegend = F)