0

s私は次のようなデータセット（pts）を持っています：

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")

次のようにして、データフレーム（pts）をマップにマップできることを知っていました。

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)

次に、さまざまな地域のptsの平均を示す別のマップを作成したいと思います。使用したいシェープファイルはftp://ftp.epa.gov/wed/ecoregions/cec_na/NA_CEC_Eco_Level2.zipからのものですが、これは北米の地図です。この北米の地図に基づいて、米国のみを表示する地図を作成するにはどうすればよいですか？または、これを行う別のより良い方法はありますか？本当にありがとう。

4

1 に答える 1

0

シェープファイルのデータだけに基づいて米国以外のデータを切り出すのは難しいと思います。これは、地域が政治的な境界に対応していないためrgeosです。

「eco」がまたはSpatialPolygonsDataFrameによって読み込まれると仮定すると、インデックス作成に使用できるキーデータを参照してください。rgdal::readOGRmaptools::readShapeSpatial

sapply(as.data.frame(eco), function(x) if(!is.numeric(x)) unique(x) else NULL)

プロットするだけの場合は、最初に米国の地域のみを含むマップを設定してから、オーバープロットします。

library(maps)
map("usa", col = "transparent")

データがランベルト正積方位図法にあることがわかります。

proj4string(eco)
[1] " +proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"

したがって、require（rgdal）eco.laea <-spTransform（eco、CRS（ "+ proj = longlat + ellpse = WGS84"））plot（eco.laea、add = TRUE）

元のランベルト正積方位図法でプロットする場合は、その投影法でバウンディングボックスを取得し、それに基づいてプロットを開始する必要があります。簡単な例として、既存のデータを使用しました。データを別の境界に対してrgeoでトリミングすることもできると確信していますが、実際に必要なものによって異なります。

于 2012-04-04T22:12:39.770 に答える