Skip to content
This repository has been archived by the owner on Jan 13, 2021. It is now read-only.

Commit

Permalink
update raster
Browse files Browse the repository at this point in the history
  • Loading branch information
Adam M. Wilson committed Oct 17, 2016
1 parent d25a024 commit ec36171
Show file tree
Hide file tree
Showing 21 changed files with 463 additions and 127 deletions.
84 changes: 69 additions & 15 deletions 05_Raster.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -509,9 +509,16 @@ ggplot(d2,aes(x=value))+
## ------------------------------------------------------------------------
transect = SpatialLinesDataFrame(
SpatialLines(list(Lines(list(Line(
cbind(c(19, 26),c(-33.5, -33.5)))), ID = "ZAF"))),
rbind(c(19, -33.5),c(26, -33.5)))), ID = "ZAF"))),
data.frame(Z = c("transect"), row.names = c("ZAF")))

# OR

transect=SpatialLinesDataFrame(
readWKT("LINESTRING(19 -33.5,26 -33.5)"),
data.frame(Z = c("transect")))


gplot(r1)+geom_tile(aes(fill=value))+
geom_line(aes(x=long,y=lat),data=fortify(transect),col="red")

Expand All @@ -524,14 +531,25 @@ gplot(r1)+geom_tile(aes(fill=value))+
trans=raster::extract(x=clim[[12:14]],
y=transect,
along=T,
cellnumbers=T)%>%
as.data.frame()
cellnumbers=T)%>%
data.frame()
head(trans)

#'
#' #### Add other metadata and reshape
## ------------------------------------------------------------------------
trans[,c("lon","lat")]=coordinates(clim)[trans$cell]
trans$order=as.integer(rownames(trans))

head(trans)

#'
## ------------------------------------------------------------------------
transl=group_by(trans,lon,lat)%>%
gather(variable, value, -lon, -lat, -cell, -order)
head(transl)

#'
## ------------------------------------------------------------------------
ggplot(transl,aes(x=lon,y=value,
colour=variable,
group=variable,
Expand Down Expand Up @@ -567,7 +585,7 @@ ggplot(rsp@data, aes(map_id = id, fill=bio1)) +


#'
#'
#' > For more details about plotting spatialPolygons, see [here](https://github.com/hadley/ggplot2/wiki/plotting-polygon-shapefiles)
#'
#' ## Example Workflow
#'
Expand All @@ -583,12 +601,32 @@ ggplot(rsp@data, aes(map_id = id, fill=bio1)) +
country=getData('GADM', country='TUN', level=1)
tmax=getData('worldclim', var='tmax', res=10)
gain(tmax)=0.1
names(tmax)

#'
#' Default layer names can be problematic/undesirable.
## ------------------------------------------------------------------------
sort(names(tmax))

## Options
month.name
month.abb
sprintf("%02d",1:12)
sprintf("%04d",1:12)

#' See `?sprintf` for details
#'
## ------------------------------------------------------------------------
names(tmax)=sprintf("%02d",1:12)

tmax_crop=crop(tmax,country)
tmaxave_crop=mean(tmax_crop) # calculate mean annual maximum temperature
tmaxavefocal_crop=focal(tmaxave_crop,
fun=median,
w=matrix(1,11,11))

#'
#' > Only a few datasets are available usig `getData()` in the raster package, but you can download almost any file on the web with `file.download()`.
#'
#' Report quantiles for each layer in a raster* object
## ------------------------------------------------------------------------
Expand All @@ -598,10 +636,9 @@ cellStats(tmax_crop,"quantile")
#'
#' ## Create a Transect (SpatialLinesDataFrame)
## ------------------------------------------------------------------------
transect = SpatialLinesDataFrame(
SpatialLines(list(Lines(list(Line(
cbind(c(8, 10),c(36, 36)))), ID = "T1"))),
data.frame(Z = c("transect"), row.names = c("T1")))
transect=SpatialLinesDataFrame(
readWKT("LINESTRING(8 36,10 36)"),
data.frame(Z = c("T1")))

#'
#'
Expand All @@ -611,16 +648,17 @@ gplot(tmax_crop)+
geom_tile(aes(fill=value))+
scale_fill_gradientn(
colours=c("brown","red","yellow","darkgreen","green"),
trans="log10",
name="Temp")+
facet_wrap(~variable)+
## now add country overlays
geom_path(data=fortify(country),
mapping=aes(x=long,y=lat,
group=group,
order=order))+
# now add transect line
geom_line(aes(x=long,y=lat),
data=fortify(transect),col="red")+
coord_equal()
data=fortify(transect),col="red",size=3)+
coord_map()

#'
#'
Expand All @@ -635,10 +673,15 @@ trans[,c("lon","lat")]=coordinates(tmax_crop)[trans$cell]
trans$order=as.integer(rownames(trans))
head(trans)


#'
#' Reformat to 'long' format.
## ------------------------------------------------------------------------
transl=group_by(trans,lon,lat)%>%
gather(variable, value, -lon, -lat, -cell, -order)%>%
separate(variable,into = c("temp","month"),4)%>%
mutate(month=as.numeric(month))
separate(variable,into = c("X","month"),1)%>%
mutate(month=as.numeric(month),monthname=factor(month.name[month],ordered=T,levels=month.name))
head(transl)

#'
#' ## Plot the transect data
Expand All @@ -654,6 +697,18 @@ ggplot(transl,
name="Month")+
geom_line()

#'
#' Or the same data in a levelplot:
## ------------------------------------------------------------------------
ggplot(transl,
aes(x=lon,y=monthname,
fill=value))+
ylab("Month")+
scale_fill_distiller(
palette="PuBuGn",
name="Tmax")+
geom_raster()

#'
#'
#' ## Raster Processing
Expand All @@ -664,4 +719,3 @@ ggplot(transl,
#' * Disk space and temporary files
#' * Use of external programs (e.g. GDAL)
#' * Use of external GIS viewer (e.g. QGIS)
#'
86 changes: 70 additions & 16 deletions 05_Raster.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -541,9 +541,16 @@ Extract values along a transect.
```{r}
transect = SpatialLinesDataFrame(
SpatialLines(list(Lines(list(Line(
cbind(c(19, 26),c(-33.5, -33.5)))), ID = "ZAF"))),
rbind(c(19, -33.5),c(26, -33.5)))), ID = "ZAF"))),
data.frame(Z = c("transect"), row.names = c("ZAF")))
# OR
transect=SpatialLinesDataFrame(
readWKT("LINESTRING(19 -33.5,26 -33.5)"),
data.frame(Z = c("transect")))
gplot(r1)+geom_tile(aes(fill=value))+
geom_line(aes(x=long,y=lat),data=fortify(transect),col="red")
```
Expand All @@ -556,14 +563,25 @@ gplot(r1)+geom_tile(aes(fill=value))+
trans=raster::extract(x=clim[[12:14]],
y=transect,
along=T,
cellnumbers=T)%>%
as.data.frame()
cellnumbers=T)%>%
data.frame()
head(trans)
```

#### Add other metadata and reshape
```{r}
trans[,c("lon","lat")]=coordinates(clim)[trans$cell]
trans$order=as.integer(rownames(trans))
head(trans)
```

```{r}
transl=group_by(trans,lon,lat)%>%
gather(variable, value, -lon, -lat, -cell, -order)
head(transl)
```

```{r}
ggplot(transl,aes(x=lon,y=value,
colour=variable,
group=variable,
Expand Down Expand Up @@ -599,7 +617,7 @@ ggplot(rsp@data, aes(map_id = id, fill=bio1)) +
```


> For more details about plotting spatialPolygons, see [here](https://github.com/hadley/ggplot2/wiki/plotting-polygon-shapefiles)
## Example Workflow

Expand All @@ -615,13 +633,33 @@ ggplot(rsp@data, aes(map_id = id, fill=bio1)) +
country=getData('GADM', country='TUN', level=1)
tmax=getData('worldclim', var='tmax', res=10)
gain(tmax)=0.1
names(tmax)
```

Default layer names can be problematic/undesirable.
```{r}
sort(names(tmax))
## Options
month.name
month.abb
sprintf("%02d",1:12)
sprintf("%04d",1:12)
```
See `?sprintf` for details

```{r}
names(tmax)=sprintf("%02d",1:12)
tmax_crop=crop(tmax,country)
tmaxave_crop=mean(tmax_crop) # calculate mean annual maximum temperature
tmaxavefocal_crop=focal(tmaxave_crop,
fun=median,
w=matrix(1,11,11))
```

> Only a few datasets are available usig `getData()` in the raster package, but you can download almost any file on the web with `file.download()`.
Report quantiles for each layer in a raster* object
```{r}
cellStats(tmax_crop,"quantile")
Expand All @@ -630,10 +668,9 @@ cellStats(tmax_crop,"quantile")

## Create a Transect (SpatialLinesDataFrame)
```{r}
transect = SpatialLinesDataFrame(
SpatialLines(list(Lines(list(Line(
cbind(c(8, 10),c(36, 36)))), ID = "T1"))),
data.frame(Z = c("transect"), row.names = c("T1")))
transect=SpatialLinesDataFrame(
readWKT("LINESTRING(8 36,10 36)"),
data.frame(Z = c("T1")))
```


Expand All @@ -643,16 +680,17 @@ gplot(tmax_crop)+
geom_tile(aes(fill=value))+
scale_fill_gradientn(
colours=c("brown","red","yellow","darkgreen","green"),
trans="log10",
name="Temp")+
facet_wrap(~variable)+
## now add country overlays
geom_path(data=fortify(country),
mapping=aes(x=long,y=lat,
group=group,
order=order))+
# now add transect line
geom_line(aes(x=long,y=lat),
data=fortify(transect),col="red")+
coord_equal()
data=fortify(transect),col="red",size=3)+
coord_map()
```


Expand All @@ -667,10 +705,15 @@ trans[,c("lon","lat")]=coordinates(tmax_crop)[trans$cell]
trans$order=as.integer(rownames(trans))
head(trans)
```

Reformat to 'long' format.
```{r}
transl=group_by(trans,lon,lat)%>%
gather(variable, value, -lon, -lat, -cell, -order)%>%
separate(variable,into = c("temp","month"),4)%>%
mutate(month=as.numeric(month))
separate(variable,into = c("X","month"),1)%>%
mutate(month=as.numeric(month),monthname=factor(month.name[month],ordered=T,levels=month.name))
head(transl)
```

## Plot the transect data
Expand All @@ -687,6 +730,18 @@ ggplot(transl,
geom_line()
```

Or the same data in a levelplot:
```{r}
ggplot(transl,
aes(x=lon,y=monthname,
fill=value))+
ylab("Month")+
scale_fill_distiller(
palette="PuBuGn",
name="Tmax")+
geom_raster()
```


## Raster Processing

Expand All @@ -695,5 +750,4 @@ Things to consider:
* RAM limitations
* Disk space and temporary files
* Use of external programs (e.g. GDAL)
* Use of external GIS viewer (e.g. QGIS)

* Use of external GIS viewer (e.g. QGIS)
Loading

0 comments on commit ec36171

Please sign in to comment.