| 1 |
e1df6e22
|
Adam M. Wilson
|
### Script to explore selection of Case Study Regions
|
| 2 |
|
|
|
| 3 |
|
|
library(rgdal)
|
| 4 |
|
|
library(latticeExtra)
|
| 5 |
|
|
library(maptools)
|
| 6 |
|
|
|
| 7 |
|
|
## make local copy of ghcn database on atlas
|
| 8 |
|
|
# pg_dump ghcn > ghcn
|
| 9 |
|
|
# copy file to local server
|
| 10 |
|
|
# createdb ghcn
|
| 11 |
|
|
# psql -d ghcn < ghcn
|
| 12 |
|
|
|
| 13 |
|
|
|
| 14 |
|
|
### get modis tiles
|
| 15 |
|
|
modt=readOGR("/home/adamw/acrobates/Global/modis_sinusoidal","modis_sinusoidal_grid_world",)
|
| 16 |
|
|
tiles=c("H18V1","H18V2","H18V3","H11V8","H19V12","H20V12","H21V9","H21V8","H31V11","H31V10","H29V5","H28V4","H28V5")
|
| 17 |
|
|
modt$roi=as.factor(modt$HV%in%tiles)
|
| 18 |
|
|
|
| 19 |
|
|
## get coastline data
|
| 20 |
|
|
#coast=readOGR("/media/data/globallayers/GSHHS/GSHHS_shp/c/","GSHHS_c_L2")
|
| 21 |
|
|
coast=Rgshhs("/media/data/globallayers/GSHHS/gshhs/gshhs_c.b", ylim=c(-60,90),xlim=c(0.01,359.99), level = 4, minarea = 1, shift = T, verbose = TRUE, checkPolygons=T)[["SP"]]
|
| 22 |
|
|
coast2=nowrapSpatialPolygons(coast)
|
| 23 |
|
|
coast=getRgshhsMap("/media/data/globallayers/GSHHS/gshhs/gshhs_c.b", ylim=c(-60,90),xlim=c(-180,180), level=4,shift = T, verbose = TRUE, checkPolygons=T)
|
| 24 |
|
|
|
| 25 |
|
|
## create union version
|
| 26 |
|
|
coast2=gUnionCascaded(coast) #dissolve any subparts of coast
|
| 27 |
|
|
coast2=SpatialPolygonsDataFrame(coast2,data=data.frame(id=1))
|
| 28 |
|
|
## create spatial lines version
|
| 29 |
|
|
coastl=as(coast,"SpatialLines")
|
| 30 |
|
|
|
| 31 |
|
|
|
| 32 |
|
|
### transform to dproj
|
| 33 |
|
|
coast=spTransform(coast,projection(modt))
|
| 34 |
|
|
coastl=spTransform(coastl,CRS(proj4string(modt)))
|
| 35 |
|
|
|
| 36 |
|
|
|
| 37 |
|
|
spplot(modt,zcol="roi",col.regions=c("transparent","red"))+layer(sp.lines(coastl))
|
| 38 |
|
|
|
| 39 |
|
|
## Bounding box of region in lat/lon
|
| 40 |
|
|
roi_ll=spTransform(roi,CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
|
| 41 |
|
|
roi_bb=bbox(roi_ll)
|
