|
1
|
### 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)
|
