サイトマップ サイト内検索 まとめ

huagrahuma.dem (topmodel パッケージ) データの紹介

このページでは，次のことを行います

前準備

使用するソフトウェア

huagrahuma, huagrahuma.dem データセット

image 関数による表示 (例) 

require(topmodel)

data(huagrahuma)
data(huagrahuma.dem)

image(huagrahuma.dem)

contour 関数による表示 (例) 

require(topmodel)

data(huagrahuma)
data(huagrahuma.dem)

contour(huagrahuma.dem)

ggplot 関数を用いた表示 (例)

【関連する外部ページ】 http://zoas.blog112.fc2.com/blog-entry-1.html 

require(ggplot2)
require(topmodel)
require(data.table)

data(huagrahuma)
data(huagrahuma.dem)


plotdem <- function(DEM, c) { 
  x <- rep(1:nrow(DEM), ncol(DEM))      # 1 2 3 1 2 3 1 2 3 1 2 3
  y <- rep(1:ncol(DEM), each=nrow(DEM)) # 1 1 1 1 2 2 2 2 3 3 3 3

  T <- data.table(x=x, y=y, val=as.numeric(DEM))
  ggplot(T, aes(x=x, y=y, fill=val)) + geom_tile() + scale_fill_gradientn("frequency", colours = c)
}

plotdem(huagrahuma.dem, cm.colors(20)) 
plotdem(huagrahuma.dem, topo.colors(20)) 
plotdem(huagrahuma.dem, terrain.colors(20)) 
plotdem(huagrahuma.dem, heat.colors(20)) 
plotdem(huagrahuma.dem, rainbow(20))

◆ plotdem(huagrahuma.dem, cm.colors(20))

◆ plotdem(huagrahuma.dem, topo.colors(20))

◆ plotdem(huagrahuma.dem, terrain.colors(20))

◆ plotdem(huagrahuma.dem, heat.colors(20))

◆ plotdem(huagrahuma.dem, rainbow(20))

persp 関数による表示 (例) 

require(topmodel)

data(huagrahuma)
data(huagrahuma.dem)

persp(huagrahuma.dem, theta = 135, phi = 30, shade=0.75, border=NA, box=FALSE )

scatterplot3d 関数による表示 (例) 

require(scatterplot3d)
require(topmodel)

data(huagrahuma)
data(huagrahuma.dem)


rgldem <- function(DEM, col) { 
  x <- rep(1:nrow(DEM), ncol(DEM))      # 1 2 3 1 2 3 1 2 3 1 2 3
  y <- rep(1:ncol(DEM), each=nrow(DEM)) # 1 1 1 1 2 2 2 2 3 3 3 3
  z <- as.numeric(DEM)
  scatterplot3d( x, y, z )
}

rgldem(huagrahuma.dem)

勾配強度

insol パッケージの cgrad + slope 

require(ggplot2)
require(topmodel)
require(insol)
require(data.table)

data(huagrahuma)
data(huagrahuma.dem)


plotdem <- function(DEM, c) { 
  x <- rep(1:nrow(DEM), ncol(DEM))      # 1 2 3 1 2 3 1 2 3 1 2 3
  y <- rep(1:ncol(DEM), each=nrow(DEM)) # 1 1 1 1 2 2 2 2 3 3 3 3

  T <- data.table(x=x, y=y, val=as.numeric(DEM))
  ggplot(T, aes(x=x, y=y, fill=val)) + geom_tile() + scale_fill_gradientn("frequency", colours = c)
}


slopedem <- function(DEM, dl_long, dl_lat) { 
  return( slope( cgrad(DEM, dl_long, dl_lat) ) )
}

plotdem(slopedem(huagrahuma.dem, 1, 1), rev(heat.colors(20)))
contour(slopedem(huagrahuma.dem, 1, 1))

◆ plotdem(slopedem(huagrahuma.dem, 1, 1), rev(heat.colors(20)))

◆ contour(slopedem(huagrahuma.dem, 1, 1))

斜面の向き

insol パッケージの cgrad + aspect 

require(ggplot2)
require(topmodel)
require(insol)
require(data.table)

data(huagrahuma)
data(huagrahuma.dem)


plotdem <- function(DEM, c) { 
  x <- rep(1:nrow(DEM), ncol(DEM))      # 1 2 3 1 2 3 1 2 3 1 2 3
  y <- rep(1:ncol(DEM), each=nrow(DEM)) # 1 1 1 1 2 2 2 2 3 3 3 3

  T <- data.table(x=x, y=y, val=as.numeric(DEM))
  ggplot(T, aes(x=x, y=y, fill=val)) + geom_tile() + scale_fill_gradientn("frequency", colours = c)
}


aspectdem <- function(DEM, dl_long, dl_lat) { 
  return( aspect( cgrad(DEM, dl_long, dl_lat) ) )
}

plotdem(aspectdem(huagrahuma.dem, 1, 1), cm.colors(20)) 
contour(aspectdem(huagrahuma.dem, 1, 1))

◆ plotdem(aspectdem(huagrahuma.dem, 1, 1), cm.colors(20))

◆ contour(aspectdem(huagrahuma.dem, 1, 1))

topmodel パッケージの sinkfill 

require(ggplot2)
require(topmodel)
require(insol)
require(data.table)

data(huagrahuma)
data(huagrahuma.dem)


plotdem <- function(DEM, c) { 
  x <- rep(1:nrow(DEM), ncol(DEM))      # 1 2 3 1 2 3 1 2 3 1 2 3
  y <- rep(1:ncol(DEM), each=nrow(DEM)) # 1 1 1 1 2 2 2 2 3 3 3 3

  T <- data.table(x=x, y=y, val=as.numeric(DEM))
  ggplot(T, aes(x=x, y=y, fill=val)) + geom_tile() + scale_fill_gradientn("frequency", colours = c)
}

plotdem( sinkfill(huagrahuma.dem, 25, 5), cm.colors(20) )
contour(sinkfill(huagrahuma.dem, 25, 5))
plotdem( sinkfill(huagrahuma.dem, 25, 8), cm.colors(20) )
contour(sinkfill(huagrahuma.dem, 25, 8))

◆ plotdem( sinkfill(huagrahuma.dem, 25, 5), cm.colors(20) )

◆ contour(sinkfill(huagrahuma.dem, 25, 5))

◆ plotdem( sinkfill(huagrahuma.dem, 25, 8), cm.colors(20) )

◆ contour(sinkfill(huagrahuma.dem, 25, 8))

topmodel パッケージの topidx 

require(ggplot2)
require(topmodel)
require(insol)
require(data.table)

data(huagrahuma)
data(huagrahuma.dem)


plotdem <- function(DEM, c) { 
  x <- rep(1:nrow(DEM), ncol(DEM))      # 1 2 3 1 2 3 1 2 3 1 2 3
  y <- rep(1:ncol(DEM), each=nrow(DEM)) # 1 1 1 1 2 2 2 2 3 3 3 3

  T <- data.table(x=x, y=y, val=as.numeric(DEM))
  ggplot(T, aes(x=x, y=y, fill=val)) + geom_tile() + scale_fill_gradientn("frequency", colours = c)
}

plotdem( topidx(huagrahuma.dem, resolution= 25)$atb, rev(heat.colors(20)) )
contour( topidx(huagrahuma.dem, resolution= 25)$atb )

◆ plotdem( topidx(huagrahuma.dem, resolution= 25)$atb, rev(heat.colors(20)) )

◆ contour( topidx(huagrahuma.dem, resolution= 25)$atb )

日照

doshade() 関数で影を求める 

require(ggplot2)
require(topmodel)
require(insol)
require(data.table)

data(huagrahuma)
data(huagrahuma.dem)


plotdem <- function(DEM, c) { 
  x <- rep(1:nrow(DEM), ncol(DEM))      # 1 2 3 1 2 3 1 2 3 1 2 3
  y <- rep(1:ncol(DEM), each=nrow(DEM)) # 1 1 1 1 2 2 2 2 3 3 3 3

  T <- data.table(x=x, y=y, val=as.numeric(DEM))
  ggplot(T, aes(x=x, y=y, fill=val)) + geom_tile() + scale_fill_gradientn("frequency", colours = c)
}


cellsize=header[5,2]
sv=normalvector(65,315)
sh=doshade(huagrahuma.dem, sv, 1)

plotdem(sh, grey(1:100/100))

関数化の途中 

## add intensity of illumination in this case sun at NW 45 degrees elevation
sv=normalvector(45, 315)
grd=cgrad(huagrahuma.dem, cellsize)
hsh=grd[,,1]*sv[1]+grd[,,2]*sv[2]+grd[,,3]*sv[3]
## remove negative incidence angles (self shading)
hsh=(hsh+abs(hsh))/2
sh=doshade(huagrahuma.dem, sv, cellsize)
plotdem(hsh*sh, grey(1:100/100))

insol パッケージに日陰、日当たりの機能あり

M <- array( c( x, y, D ), dim=c(nrow(D), ncol(D), 3) );

sunpos() Examples ## Julian Day hourly intervals at spring equinox jd=JD(seq(ISOdate(2012,3,20,0),ISOdate(2012,3,20,23),by="hour")) ## sun position sp=sunpos(sunvector(jd,46.813,9.844,1)) ## daylight zenith<=90 sp=sp[which(sp[,2]<=90),] ## Plot the apparent solar path at Davos on the spring equinox ramp = colorRamp(c("red", "orange","yellow")) crmp=c(rgb(ramp(seq(1/6,1,1/6)), max = 255),rgb(ramp(seq(1,1/6,-1/6)), max = 255)) plot(sp[,1],90-sp[,2],xlab='Azimuth', ylab='Elevation',main='Apparent solar path at Davos on the spring equinox', pch=20,col=crmp,cex=(300-sp[,2])/90) ## Not run: require(plotrix) polar.plot(90-sp[,2],sp[,1],start=90,clockwise=TRUE,rp.type='s', point.symbols=20,point.col=2,cex=2,radial.lim=c(0,90), main='Apparent solar path at Davos on the spring equinox') ## End(Not run) sunvector Solar vector Description Calculates a unit vector in the direction of the sun from the observer position. Usage sunvector(jd, latitude, longitude, timezone)sunvector Arguments jd Julian Day and decimal fraction. latitude Latitude of observer in degrees and decimal fraction. longitude Longitude of observer in degrees and decimal fraction. timezone Time zone, west is negative. Details To calculate the sunvector to the nearest hour, give the Julian Day with a precission better than 1/24; to approximate it to the nearest minute use decimal fractions smaller than 1/(24*60), and so on. Value 3 column matrix with the x, y , z coordinates of the sun vector. Author(s) Javier G. Corripio References Corripio, J. G.: 2003, Vectorial algebra algorithms for calculating terrain parameters from DEMs and the position of the sun for solar radiation modelling in mountainous terrain, International Journal of Geographical Information Science 17(1), 1-23. sunpos Examples # Current solar vector at Greenwich observatory sunvector(JD(Sys.time()),51.4778,-0.0017,0) juneday=JD(seq(ISOdate(2012,6,21,0),ISOdate(2012,6,21,23,30),by=’30 min’)) ## Not run: # Path of the sun over Greenwich in summer require(scatterplot3d) scatterplot3d(sunvector(juneday,51.4778,-0.0017,0), ylim=c(-1,1),zlim=c(0,1),pch=8,color=’orange’) ## End(Not run) # print values options(digits=12) # make sure decimals are printed sunvector(juneday,51.4778,-0.0017,0)34

以下、書きかけ

topmodel パッケージの topmodel 

data(huagrahuma)
attach(huagrahuma)
## returns the simulated runoff (Qobs not given)
Qsim <- topmodel(parameters, topidx, delay, rain,ET0, verbose = TRUE)
## plot observed and simulated discharge:
plot(Qobs)
points(Qsim$Q, col="red", type="l")