2012-03-29 4 views
1

Ich habe den Datensatz (pts) wie folgt aus:Wie Polygon zu berechnen bedeutet und sie zuordnen?

x <- seq(-124.25,length=115,by=0.5)  
y <- seq(26.25,length=46,by=0.5) 
z = 1:5290 

longlat <- expand.grid(x = x, y = y) # Create an X,Y grid 
pts=data.frame(longlat,z) 
names(pts) <- c("x","y","data") 

Ich wusste, dass ich die Datenrahmen (pts) in eine Karte, indem Sie abbilden:

library(sp) 
library(rgdal) 
library(raster) 
library(maps) 
coordinates(pts)=~x+y 

proj4string(pts)=CRS("+init=epsg:4326") # set it to long, lat 

pts = spTransform(pts,CRS(" +init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")) 
pts <- as(pts, "SpatialPixelsDataFrame") 
r = raster(pts) 
projection(r) = CRS(" +init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0") 

plot(r) 
map("usa",add=T) 

meine Frage ist, kann ich das, wie berechnen die Mittel der 10 EPA-Regionen und die Mittel Karte?

die WPA-Regionen am unteren Ende der Webseite kann http://www.epa.gov/wed/pages/ecoregions/level_iii_iv.htm

bei finden

Dank

+0

Plotten Bevor auch auf Ihre Frage bekommen, gibt es mehrere Probleme mit Ihrem Code. Ich habe einen behoben, indem ich einen Aufruf zu 'library (raster)' hinzugefügt habe. Wenn Sie möchten, dass "r" die Werte in "z" enthält, müssen Sie vor dem Aufruf von "raster (pts)" auch die folgende Zeile hinzufügen: pts <- as (pts, "SpatialPixelsDataFrame"). (Weil Sie dies nicht getan haben, schlägt 'plot (r)' zur Zeit mit einer Fehlermeldung fehl). Um dies zu einem leicht reproduzierbaren Beispiel zu machen, können Sie auch Code hinzufügen, um ein einfaches SpatialPolygon-Objekt zu erstellen, das den RasterLayer überlagert. (Sonst müssen die Leute einen erstellen, um Ihnen zu helfen). –

+0

Ich kann rgdal nicht mit R 2.14.2 arbeiten. Es gibt ein Problem mit dem Namespace, das entweder anzeigt, dass ich ein Problem habe oder das Paket nicht aktualisiert wurde und Sie eine ältere Version von R verwenden. Die beste Hilfe, die ich geben kann, ist ein Link zu einer Frage, die ich gestellt habe das kann ähnlich sein [(LINK)] (http://stackoverflow.com/questions/9441778/improve-centering-county-names-ggplot-maps) –

+0

@TylerRinker - Ich renne auch R 2.14.2, auf Windows XP, sowie die neueste Version von ** rgdal ** ('0.7.8', veröffentlicht 2012-01-18), und es funktioniert ohne Probleme. Aus Neugier, sind Sie auf Windows oder etwas anderes? –

Antwort

2

zuerst die Formdatei lesen

er <- readOGR("Eco_Level_III_US.shp", "Eco_Level_III_US") 

Dann stellen Sie sicher, ther Raster r und die Ökoregionen er die gleiche Projektion

er.4326 <- spTransform(er, CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")) 

die r Rasterdaten aus dem Shape-Datei extrahieren (das kann ein paar Minuten dauern), berechnen Sie den Mittelwert und fügen Sie ihn zu Ihren Polyongs hinzu.

er.v <- extract(r, er.4326) 
means <- sapply(er.v, mean) 
er.4326$means <- means 

Und es schließlich

spplot(er.4326, "means") 
+0

Ehrfürchtig. Wie wäre es möglich, für jede Region eine Legende zu erstellen, die der auf ihrer ursprünglichen Karte gezeigten entspricht? Entschuldige die lahmen Fragen. – Dan

+0

Ich habe festgestellt, dass es 31 Regionen auf dieser Karte gibt. Gibt es eine Möglichkeit, die Mittelwerte der Polygone zu berechnen, die auf der Karte unten auf der Webseite angezeigt werden? Ich konnte für diese Karte kein Shapefile finden. – Dan

Verwandte Themen