#
# in UNIX untar all the tiles:
# e.g. 
# tar -z -xvf w140n40.tar.gz ; rm *.DMW *.PRJ *.SCH *.SRC *.STX


library( fields)
source("get.GTOPO30.R")

#> range( lon)-360
#[1] -147.86960  -44.30399
#> range( lat)
#[1] 21.84493 67.17414

xr<- c( -148-15, -44+15)
yr<- c( 21- 5, 67 + 10)

# check
#plot( xr, yr, type="n"); world( add=TRUE)
#points( lon-360, lat, pch=".", col="red")

tiles<- c(
"W180N40",
"W140N40",
"W100N40",
"W060N40",

"W180N90",
"W140N90",
"W100N90",
"W060N90")

get.GTOPO30( tiles[1], shrink =TRUE)-> tile1
get.GTOPO30( tiles[2], shrink =TRUE)-> tile2
get.GTOPO30( tiles[3], shrink =TRUE)-> tile3
get.GTOPO30( tiles[4], shrink =TRUE)-> tile4
get.GTOPO30( tiles[5], shrink =TRUE)-> tile5
get.GTOPO30( tiles[6], shrink =TRUE)-> tile6
get.GTOPO30( tiles[7], shrink =TRUE)-> tile7
get.GTOPO30( tiles[8], shrink =TRUE)-> tile8

plot( xr,yr, type="n")
image( tile1, add=TRUE)
image( tile2, add=TRUE)
image( tile3, add=TRUE)
image( tile4, add=TRUE)
image( tile5, add=TRUE)
image( tile6, add=TRUE)
image( tile7, add=TRUE)
image( tile8, add=TRUE)

# check alignment

 table(round(tile1$x-tile2$x,7))
 table(round(tile1$x-tile3$x,7))
 table(round(tile1$x-tile4$x,7))
 table(round(tile1$x-tile5$x,7))
 table(round(tile1$x-tile6$x,7))
 table(round(tile1$x-tile7$x,7))

 table(round(tile1$y-tile2$y,7))
 table(round(tile1$y-tile3$y,7))
 table(round(tile1$y-tile4$y,7))
 table(round(tile1$y-tile5$y,7))
 table(round(tile1$y-tile6$y,7))
 table(round(tile1$y-tile7$y,7))

# check image

test<- list( x= c(tile1$x, tile2$x, tile3$x, tile4$x), 
         y= tile2$y, z= rbind( tile1$z,tile2$z, tile3$z, tile4$z) ) 

test2<- list( x= c(tile5$x, tile6$x, tile7$x, tile8$x), 
         y= tile5$y, z= rbind( tile5$z,tile6$z, tile7$z, tile8$z) ) 


NARCCAP.elev<- list( x= test$x, y= c(test$y, test2$y), 
                       z= cbind(test$z, test2$z) )
 
save( NARCCAP.elev, file="NARCCAP.elev.rda")

# testing image
#
load("NARCCAP.elev.rda")

png( "NARCCAP.elev.png", height=320, width=320)
image.plot(NARCCAP.elev, xlab="", ylab="", col=terrain.colors(256))
dev.off()

png( "NARCCAP.elev.high.png", height=4*320, width=4*320)
image.plot(NARCCAP.elev, xlab="", ylab="", col=terrain.colors(256))
dev.off()