---
title: day 1 practical session, spatial data in R
---

If your map is upside-down, you'll need the new version of the mapmisc package
```{r installMapmics, eval=FALSE}
install.packages("mapmisc", 
	repos="http://r-forge.r-project.org")
```


```{r packages}
library('mapmisc')
```

# A nice map

consider the following code.

```{r patrick, fig.height=4, fig.width=5}
mycoords = rbind(c(-2.785556, 54.010278), c(-2.736, 53.761 ))
patrick = vect(mycoords, 
	atts=data.frame(name=c('Work','Home')),
	crs='epsg:4326') # long-lat
#patrick = project(patrick, 
#	mapmisc::tpers(patrick[2:1,], 
#		hKm = 600*1000))	
patrick = project(patrick, "epsg:27700") # uk national grid
theMap = openmap(patrick, 
	buffer=1000*1000*c(15,15,5,5), fact=2) # 15km extra left,right, 5km north, south
map.new(theMap)
plot(theMap, add=TRUE)
plot(patrick, add=TRUE, col='blue')
mapmisc::scaleBar(crs(patrick), 'topleft', bg='white', cex=2)
forInset = openmap(patrick, buffer=500*1000, crs=mapmisc::crsMerc)
mapmisc::insetMap(crs(patrick), 'topright', map=forInset, 
	inset=0, width=0.4)
```





It makes a map showing where I lived during my PhD years, and the location of Lancaster University.

- Choose two locations important to you, such as the location of your first primary school and where you lived at the time. 
- Find their coordinates, using google maps or openstreetmap.org
- make a map similar to mine.
- use the `text` function to label the points



# Make a map involving polygons

Use your dataset of interest to make a map involving polygons.  Or if you prefer do the following.

On the web site [sites.wustl.edu/acag/datasets/surface-pm2-5/#V6.GL.02.04](https://sites.wustl.edu/acag/datasets/surface-pm2-5/#V6.GL.02.04) there are images of particulate matter concentrations worldwide.  I downloaded the 2020 global image from [wustl.app.box.com/s/y143mciw7jz7ft2qe3hccjw65m3xe8f2/folder/327753085334](https://wustl.app.box.com/s/y143mciw7jz7ft2qe3hccjw65m3xe8f2/folder/327753085334), saved it in my Downloads folder and 

```{r data}
myPath = normalizePath('~/Downloads')
x= rast(file.path(myPath, 
	"V6GL02.04.CNNPM25.GL.202001-202012.nc"))
brazil = project(
	geodata::gadm("BRA", level=1, path=myPath),
	crs(x))
xCrop = crop(x, ext(brazil))
xCrop[xCrop<= 0] <- NA

myCol = mapmisc::colourScale(xCrop, breaks=9, style='equal', col='RdYlGn', rev=TRUE, transform='sqrt',dec=0)

map.new(xCrop)
plot(xCrop, col=myCol$col, breaks=myCol$breaks, legend=FALSE, add=TRUE)
plot(brazil, add=TRUE, border='red')
legendBreaks('bottomright', myCol)
```

- Find the average PM by admin level 2 in Brazil.  How does Rio rank?
- Look at the help files for geodata::world.  Use it to create a similar map of all south america.

```{r pmByMuiciaplity}
brazil2 = geodata::gadm("BRA", level=2, 
	path=myPath)
brazil3 = project(brazil, crs(x))
thePm = extract(brazil3, x)	
```


# Mapping a raster

Make a map of a raster related to your project.  Or if you prefer use one of the datasets from the following web sites measure nighttime light using satelites.  

- https://gee-community-catalog.org/projects/hntl/
- https://figshare.com/articles/dataset/Harmonization_of_DMSP_and_VIIRS_nighttime_light_data_from_1992-2018_at_the_global_scale/9828827/5

- Choose an interesting part of the planet (your hometown? North Korea? Tierra del Fuego?) and plot the lights image.  
- put one or two interesting cities or places on the map





