私は、薄板スプライン アルゴリズムを使用して英国のグリッド降雨データを生成し、R の陸地にない値を排除しようとしています。これは、これまでのところ手動でしか達成できないプロセスです。この問題は (私にとっては) 挑戦的であり、説明するのも困難です。どんな助けでも大歓迎です。
最初に、R にデータ テーブルを読み込みます。このデータ テーブルは、多数のポイント ロケーションの気象観測所からの 1 日の降雨量を表しており、データ テーブルの各行には、日付、観測所の ID、観測所の東向きと北向き、その場所の日降水量と年間の平均降水量。ライブラリ フィールド、maptools、および gstat もロードします。
library(fields)
library(maptools)
library(gstat)
dat <- read.table("1961month1day1.csv", header=T, sep=",", quote = "")
names(dat) <- c("easting", "northing", "dailyrainfall","avaerageyearlyrainfall")
データのサンプルを次に示します。
dput(head(dat, 20))
structure(list(easting = c(130000L, 145000L, 155000L, 170000L,
180000L, 180000L, 180000L, 180000L, 185000L, 200000L, 200000L,
205000L, 210000L, 220000L, 225000L, 230000L, 230000L, 230000L,
230000L, 235000L), northing = c(660000L, 30000L, 735000L, 40000L,
30000L, 45000L, 60000L, 750000L, 725000L, 50000L, 845000L, 65000L,
770000L, 105000L, 670000L, 100000L, 620000L, 680000L, 95000L,
120000L), dailyrainfall = c(9.4, 4.1, 12.4, 2.8, 1.3, 3.6, 4.8, 26.7, 19.8,
4.6, 1.7, 4.1, 12.7, 1.8, 3, 5.3, 1, 1.5, 1.5, 4.6), averageyearlyrainfall = c(1334.626923,
1123.051923, 2072.030769, 1207.584615, 928, 1089.334615, 880.0884615,
2810.323077, 1933.719231, 1215.642308, 2644.171154, 1235.913462,
2140.111538, 1010.436538, 1778.432692, 1116.934615, 912.2807692,
1579.386538, 1085.498077, 1250.601923)), .Names = c("easting",
"northing", "dailyrainfall", "averageyearlyrainfall"), row.names = c(NA, 20L), class = "data.frame")
次に、薄板スプラインをデータに適合させて、グリッド サーフェスを作成し、サーフェスをプロットします。
fit <- Tps(cbind(dat$easting,dat$northing),dat$dailyrainfall)
surface(fit)
次に、以下を使用して 1 km 刻みで英国のグリッドを作成できます。
xvals <- seq(0, 700000, by=1000)
yvals <- seq(0, 1250000, by=1000)
次に、サーフェスをこのグリッドにプロットし、データをテーブルに書き込みます。
griddf <- expand.grid(xvals, yvals)
griddf$pred <- predict(fit, x=as.matrix(griddf))
write.table(griddf, file="1Jan1961grid.csv", sep=",", qmethod="double")
素晴らしい-これまでのところとても良い. これで、ポイント データを 0 ～ 700000 (E) および 0 ～ 1250000 (N) グリッド全体の 1 km グリッド データに変換しました。書き込まれたデータ テーブルは、インデックス、東座標、北座標、および降雨量の予測値を含むリストです。
ここでの課題 - このリストから、土地を超えていない値を削除したいと思います。データを Excel (または Access) に読み込み、同じグリッドと年間平均降雨量を含む別のファイル (ファイルは 1kmgridaveragerainfall.csv と呼ばれます) とデータを比較することで、これを手動で行うことができます。このファイルのサンプルを次に示します。
dput(head(dat1, 20))
structure(list(easting = c(-200000L, -200000L, -200000L, -200000L,
-200000L, -200000L, -200000L, -200000L, -200000L, -200000L, -200000L,
-200000L, -200000L, -200000L, -200000L, -200000L, -200000L, -200000L,
-200000L, -200000L), northing = c(1245000L, 1240000L, 1235000L,
1230000L, 1225000L, 1220000L, 1215000L, 1210000L, 1205000L, 1200000L,
1195000L, 1190000L, 1185000L, 1180000L, 1175000L, 1170000L, 1165000L,
1160000L, 1155000L, 1150000L), averageyearlyrainfall = c(-9999, -9999, -9999,
-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999,
-9999, -9999, -9999, -9999, -9999, -9999, -9999, -9999)), .Names = c("easting",
"northing", "averageyearlyrainfall"), row.names = c(NA, 20L), class = "data.frame")
土地の上にないグリッド スクエアの年間平均降水量は -9999 です。したがって、一致すると (つまり、Access で vlookup またはクエリを使用して)、この -9999 値を持つ値を除外できます。これにより、土地の値のみの東向きと北向き、毎日の降雨量、および年間の平均降雨量を含むデータ テーブルが残ります。次に、これを R にロードして、次を使用してプロットできます。
quilt.plot(cbind(dat$easting,dat$northing),dat$mm, add.legend=TRUE, nx=654, ny=1209,xlim=c(0,700000),ylim=c(0,1200000))
そして、英国の陸地（海域ではなく）の降雨量のプロットが残っています。
それで、誰かが同じことを達成する方法を提案できますが、Excelまたはアクセスを使用してすべてのフィルタリングなどを行わずに、つまりRのみを使用して同じことを達成できますか? 最初に両方のデータ テーブルを R にロードし、ポイント データの TPS を平均データに合わせて、-9999 に等しいグリッド スクエアがプロットされないようにする方法はありますか。
共変量 (Z) を使用して TPS を重み付けできることは知っていますが、これはまったく役に立ちますか? すなわち
fit <- Tps(cbind(dat$easting,dat$northing),dat$dailyrainfall, Z=dat$averageyearlyrainfall)
また、元のTPSの表面（フィット）を実行するとき、プロットをプロットの端まで拡張するにはどうすればよいですか-interp = TRUEのようなものをどこかに置いた場所でこれを読んだと確信していますが、これはそうではありません仕事。
どんな助けでも大歓迎です
ありがとう、トニー