Traçage de données interpolées sur la carte

J'ai des données de l'enquête de la richesse en espèces qui a été pris à différents endroits dans la Baie de Chesapeake, etats-unis, et je tiens à présenter graphiquement les données comme une "carte de chaleur."

J'ai un dataframe de coordonnées lat/long et la richesse des valeurs, que j'ai converti en un SpatialPointsDataFrame et utilisé le autoKrige() fonction de l'automap paquet pour générer les valeurs interpolées.

Tout d'abord, quelqu'un peut-il dire si je suis correctement mise en œuvre de la autoKrige() fonction?

Deuxième, j'ai du mal à la représentation des données et en y ajoutant une carte de la région. Sinon, pourrais-je spécifier l'interpolation de la grille afin de refléter les frontières de la Baie (comme l'a suggéré ici)? Des idées sur comment je pourrais faire et où je pourrais obtenir cette information? La fourniture de la grille de autoKrige() semble assez facile.

EDIT: Merci à Paul pour son super utile post! Voici ce que j'ai maintenant. La difficulté à obtenir les ggplot à accepter à la fois les données interpolées et la projection de la carte:

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")
  • l'exemple de code donne-moi un vide parcelle avec la ggplot_0.9.3 dans R 2.15 et une erreur de la parcelle R 3.0.0 (pas de méthode applicable pour la profondeur)
InformationsquelleAutor jslefche | 2012-04-06
  1. 47

    J'ai un certain nombre de remarques sur votre poste:

    À l'aide du krigeage

    Je vois que vous êtes à l'aide de la géostatistique pour construire votre heatmap. Vous pourriez aussi envisager d'autres techniques d'interpolation comme les splines (par ex. plaque Mince splines dans les domaines de package). Ceux-ci font de moins en moins des hypothèses sur les données (par exemple stationnarité), et peut également visualiser vos données de l'amende juste. La réduction du nombre d'hypothèses qui pourraient aider dans le cas où vous l'envoyez à un journal, alors vous avez moins d'expliquer à l'ensemble des évaluateurs. Vous pouvez également comparer quelques techniques d'interpolation si vous le souhaitez, voir un rapport que j'ai écrit pour quelques conseils.

    De projection de données

    Je vois que vous êtes en utilisant la latitude et de la longitude coordonnées pour le krigeage. Edzer Pebesma (auteur de gstat) a fait remarquer que il n'y a pas variogramme des modèles qui sont adaptés pour la latitude /longitude coordonnées. C'est parce que dans lat lon les distances ne sont pas directement (c'est à dire Euclidienne), mais sur une sphère, (c'est à dire Grand cercle distances). Il n'y a pas de covariance des fonctions (ou des modèles de variogramme) qui sont valables pour les coordonnées sphériques. Je recommande de les projeter à l'aide de spTransform de la rgdal paquet avant de l'utiliser automap.

    La rgdal paquet utilise la proj4 bibliothèque de projection pour effectuer les calculs. Le projet de vos données, vous devez d'abord définir sa projection:

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

    La proj4 chaîne sur le côté droit de l'expression ci-dessus définit le type de projection (+proj), le ellips qui a été utilisé (+ellps) et la référence (+datum). Pour comprendre la signification de ces termes, vous devez imaginer la Terre comme une pomme de terre. La Terre n'est pas parfaitement sphérique, il est défini par le ellips. Ni est la Terre un parfait ellipsoïde, mais la surface est plus irrégulière. Cette irrégularité est défini par la donnée. Voir aussi cet article sur Wikipédia.

    Une fois que vous avez la projection défini, vous pouvez utiliser spTransform:

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

    où les CRS("+proj etc") définit l'objectif de projection. Qui de projection dépend de votre situation géographique et de la taille de la zone d'étude.

    Comploter avec ggplot2

    Pour ajouter des polygones ou des polylignes à ggplot, s'il vous plaît un coup d'oeil à la documentation de coord_map. Cela inclut un exemple d'utilisation de la maps paquet de tracer les frontières des pays. Si vous avez besoin de charger par exemple les fichiers de formes pour votre zone d'étude, vous pouvez le faire en utilisant rgdal. N'oubliez pas que ggplot2 fonctionne avec des données.du cadre, pas SpatialPolygons. Vous pouvez transformer SpatialPolygons à data.frame à l'aide de:

    poly_df = fortify(poly_Spatial)

    Voir aussi cette fonction j'ai créé pour tracer spatiale des grilles. Il fonctionne directement sur SpatialGrids/Pixels. Notez que vous avez besoin de source d'un ou de deux autres fichiers à partir de ce référentiel (continuousToDiscrete).

    La création d'interpolation de la grille de

    J'ai créé automap pour générer une sortie de la grille, lorsqu'il n'a été spécifié. Ceci est fait par la création d'une enveloppe convexe autour des points de données et d'échantillonnage de 5000 points à l'intérieur d'elle. Les limites de la prédiction de la région, et le nombre de points échantillonnés (et donc de la résolution) est tout à fait arbitraire. Pour une application spécifique, la forme de la prédiction de la zone peut être dérivée à partir d'un polygone à l'aide de spsample de points d'échantillonnage à l'intérieur du polygone. Combien de points de l'échantillon, et donc la résolution, dépend de deux choses:

    • le type de données que vous avez, Par exemple, si vos données est très lisse, il n'y a pas beaucoup de point dans la sensibilisation de la résolution très élevé en comparaison de cette fluidité. Sinon, si vos données sont nombreux petits strcutures, vous avez besoin d'une haute résolution. Cela n'est possible bien sûr, si vous avez des observations à l'appui de cette haute résolution.
    • la densité de données. Si vos données est de plus en plus dense, vous pouvez augmenter la résolution.

    Si vous utilisez votre carte interpolée pour les analyses ultérieures, l'obtention de la résolution est importante. Si vous utilisez la carte purement pour visuatlisation fins, c'est moins important. Notez cependant que dans les deux cas, une trop haute résolution peut être trompeuse quant à la précision de vos prévisions, et qu'une trop faible résolution ne rend pas justice aux données.

    • Nous sommes dans le cou de Paul du bois, maintenant! +1 pour un bien écrit (et obtenu) réponse.
    • Extraordinairement bien écrit et instructif réponse! Merci beaucoup pour la fourniture de tels détails pour quelqu'un (moi) qui est (évidemment) très nouveau pour cette. J'ai une question: où puis-je trouver les arguments qui seraient appropriés pour CRS() dans spTransform pour la balance que je suis en train de travailler avec? (par exemple, la Baie de Chesapeake: 300 x 50 km au nord-ouest de la Mi-région de l'Atlantique)
    • Vous pouvez akways utiliser l'UTM, qui couvre l'ensemble du globe. Il le fait en se projetant dans les zones à travers le monde. Une fois que vous connaissez le numéro de zone, vous êtes bon pour aller. Pour le proj4string, il suffit de google utm proj4.Alternativement, vous pouvez utiliser un système de coordonnées local. Googler pour vous pays ou d'une zone plus coordinate system devrait. Ou de regarder ce que les cartes existantes de la zone d'utilisation.
    • Super! Supposons que j'utilise +proj=UTM +zone=18s. Je suis maintenant de la difficulté à obtenir ggplot à accepter à la fois les données interpolées et les données de la carte. Voir édité mon post ci-dessus
    • Ggplot2 utilise une projection différente de la bibliothèque de rgdal autant que je sache. Donc, mercator n'est pas égal à l'UTM, qui est un spécial mercator. L'idée est qu'une fois que toutes vos données (polygone d'une grille) sont projetés dans la même projection à l'aide de proj4 vous n'avez pas besoin de fournir des informations supplémentaires à ggplot.
    • Votre meilleur pari est de créer la grille dans la bonne projection et de changer le polygone à cette projection. Des une grille conduit à ré-échantillonnage, qui rgdal ne pas le faire.
    • Roger qui. Bon à savoir que PROJ.4 projections ne sont pas les mêmes que celles de la coord_map. Je vais convertir en données.cadre à l'aide de rgdal et de mettre à jour le post (avec code et de la parcelle) pour que les autres la référence en bas de la ligne. Merci encore, Paul!

    InformationsquelleAutor Paul Hiemstra