24

米国チェサピーク湾のさまざまな場所で採取された種の豊富さの調査データがあり、そのデータを「ヒート マップ」としてグラフィカルに提示したいと考えています。

緯度/経度座標と豊かさの値のデータフレームがあり、これを に変換し、automap パッケージ SpatialPointsDataFrameの関数を使用して補間値を生成しました。autoKrige()

まず、autoKrige()関数を正しく実装しているかどうかについてコメントできますか?

次に、データをプロットして地域の地図を重ねるのに問題があります。または、ベイの境界を反映するように補間グリッドを指定できますか (ここで提案されているように)? どうすればそれができるか、どこでその情報を入手できるかについて何か考えはありますか? グリッドを提供するのautoKrige()は簡単に見えます。

編集: 非常に役立つ投稿をしてくれた Paul に感謝します! これが私が今持っているものです。ggplot が補間されたデータと地図投影の両方を受け入れるのに問題があります。

require(rgdal)
require(automap)
#Generate lat/long coordinates and richness data
set.seed(6)
df=data.frame(
  lat=sample(seq(36.9,39.3,by=0.01),100,rep=T),
  long=sample(seq(-76.5,-76,by=0.01),100,rep=T),
  fd=runif(10,0,10))
initial.df=df

#Convert dataframe into SpatialPointsDataFrame
coordinates(df)=~long+lat

#Project latlong coordinates onto an ellipse
proj4string(df)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
#+proj = the type of projection (lat/long)
#+ellps and +datum = the irregularity in the ellipse represented by planet earth

#Transform the projection into Euclidean distances
project_df=spTransform(df, CRS("+proj=merc +zone=18s +ellps=WGS84 +datum=WGS84")) #projInfo(type="proj")

#Perform the interpolation using kriging
kr=autoKrige(fd~1,project_df)
#Extract the output and convert to dataframe for easy plotting with ggplot2
kr.output=as.data.frame(kr$krige_output)
#Plot the output
#Load the map data for the Chesapeake Bay
cb=data.frame(map("state",xlim=range(initial.df$long),ylim=range(initial.df$lat),plot=F)[c("x","y")])

ggplot()+
  geom_tile(data=kr.output,aes(x=x1,y=x2,fill=var1.pred))+  
  geom_path(data=cb,aes(x=x,y=y))+
  coord_map(projection="mercator")
4

1 に答える 1

47

あなたの投稿にはいくつかのコメントがあります：

クリギングの使用

ヒートマップの作成に地球統計学を使用しているようです。スプラインなどの他の補間手法を検討することもできます (例: fields パッケージの薄板スプライン)。これらは、データに関する仮定 (定常性など) を少なくし、データを適切に視覚化することもできます。仮定の数を減らすことは、ジャーナルに投稿する場合に役立つ可能性があります。その場合、査読者に説明する必要が少なくなります。必要に応じて、いくつかの補間手法を比較することもできます。いくつかのヒントについては、私が書いたレポートを参照してください。

データ投影

クリギングに緯度経度座標を使用しているようです。Edzer Pebesma (の著者gstat)は、緯度経度座標に適したバリオグラム モデルはないと述べています。これは、緯度経度では距離が直線 (つまり、ユークリッド) ではなく、球上 (つまり、大円距離) であるためです。球座標に有効な共分散関数 (またはバリオグラム モデル) はありません。automap を使用する前にspTransform、パッケージから使用して投影することをお勧めします。rgdal

rgdal パッケージは、proj4 投影ライブラリを使用して計算を実行します。データを射影するには、まずその射影を定義する必要があります。

proj4string(df) = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

上記の式の右側にある proj4 文字列は、投影のタイプ ( +proj)、使用された楕円 ( +ellps)、およびデータム ( +datum) を定義します。これらの用語が何を意味するかを理解するには、地球をジャガイモとして想像する必要があります。地球は完全な球形ではなく、これは楕円によって定義されます。地球も完全な楕円体ではありませんが、表面はより不規則です。この不規則性はデータムによって定義されます。ウィキペディアのこの記事も参照してください。

投影を定義したら、次を使用できますspTransform。

project_df = spTransform(df, CRS("+proj= etcetc"))

ここで、CRS("+proj etc") はターゲット投影を定義します。どの図法が適切かは、地理的な場所と調査範囲の大きさによって異なります。

ggplot2 によるプロット

ポリゴンまたはポリラインを ggplot に追加するには、 のドキュメントを参照してくださいcoord_map。mapsこれには、パッケージを使用して国の境界をプロットする例が含まれています。たとえば、調査エリアのシェープファイルをロードする必要がある場合は、 を使用して実行できますrgdal。ggplot2ではなく、data.frame で動作することを覚えておいてくださいSpatialPolygons。SpatialPolygons以下をdata.frame使用して変換できます。

poly_df = fortify(poly_Spatial)

空間グリッドをプロットするために作成したこの関数も参照してください。SpatialGrids/Pixels で直接動作します。そのリポジトリから 1 つまたは 2 つの追加ファイルを入手する必要があることに注意してください ( continuousToDiscrete )。

補間グリッドの作成

何も指定されていないときに出力グリッドを生成する automap を作成しました。これは、データ ポイントの周囲に凸包を作成し、その中の 5000 ポイントをサンプリングすることによって行われます。予測領域の境界、およびその領域でサンプリングされるポイントの数 (したがって解像度) は、まったく恣意的です。spsample特定のアプリケーションでは、予測領域の形状は、多角形内のポイントをサンプリングするために使用して、多角形から導き出すことができます。サンプリングするポイントの数、つまり解像度は、次の 2 つの要素に依存します。

  • たとえば、データが非常に滑らかな場合、この滑らかさに比べて解像度を非常に高くしてもあまり意味がありません。または、データに小規模な構造が多数ある場合は、高解像度が必要です。これはもちろん、この高解像度をサポートする観測がある場合にのみ可能です。
  • データの密度。データの密度が高い場合は、解像度を上げることができます。

補間されたマップを後続の分析に使用する場合、適切な解像度を取得することが重要です。純粋に視覚化の目的でマップを使用する場合、これはそれほど重要ではありません。ただし、どちらの場合も、解像度が高すぎると予測の精度に関して誤解を招く可能性があり、解像度が低すぎるとデータが正当化されないことに注意してください。

于 2012-04-07T22:18:41.723 に答える