15

これを行う方法のかなり単純な例を見つけましたが、うまくいきません。私はRにかなり慣れていません

library(rgdal) 
xy <- cbind(c(118, 119), c(10, 50)) 
project(xy, "+proj=utm +zone=51 ellps=WGS84") 
          [,1]    [,2] 
[1,] -48636.65 1109577 
[2,] 213372.05 5546301

しかし、これは例の数字です。変換する必要がある何千もの座標があり、それらをテーブルからこのスクリプトに取得する方法がわかりません

私のデータセットには ID、X、Y の 3 つの列があります。この式を使用してそれらを変換するにはどうすればよいですか? 私は何週間もこれにこだわっています

4

5 に答える 5

26

適切な投影メタデータが座標に関連付けられたすべてのステップにあることを確認するには、SpatialPointsDataFrameできるだけ早くポイントをオブジェクトに変換することをお勧めします。

?"SpatialPointsDataFrame-class"単純な data.frames またはマトリックスをSpatialPointsDataFrameオブジェクトに変換する方法の詳細については、を参照してください。

library(sp)
library(rgdal)

xy <- data.frame(ID = 1:2, X = c(118, 119), Y = c(10, 50))
coordinates(xy) <- c("X", "Y")
proj4string(xy) <- CRS("+proj=longlat +datum=WGS84")  ## for example

res <- spTransform(xy, CRS("+proj=utm +zone=51 ellps=WGS84"))
res
#            coordinates ID
# 1 (-48636.65, 1109577)  1
# 2    (213372, 5546301)  2

## For a SpatialPoints object rather than a SpatialPointsDataFrame, just do: 
as(res, "SpatialPoints")
# SpatialPoints:
#              x       y
# [1,] -48636.65 1109577
# [2,] 213372.05 5546301
# Coordinate Reference System (CRS) arguments: +proj=utm +zone=51
# +ellps=WGS84
于 2013-09-05T16:01:48.103 に答える
19

緯度と経度のポイントを UTM に変換する

library(sp)
library(rgdal)

#Function
LongLatToUTM<-function(x,y,zone){
 xy <- data.frame(ID = 1:length(x), X = x, Y = y)
 coordinates(xy) <- c("X", "Y")
 proj4string(xy) <- CRS("+proj=longlat +datum=WGS84")  ## for example
 res <- spTransform(xy, CRS(paste("+proj=utm +zone=",zone," ellps=WGS84",sep='')))
 return(as.data.frame(res))
}

# Example
x<-c( -94.99729,-94.99726,-94.99457,-94.99458,-94.99729)
y<-c( 29.17112, 29.17107, 29.17273, 29.17278, 29.17112)
LongLatToUTM(x,y,15)
于 2015-05-13T21:52:21.420 に答える
6

あなたの質問では、データセットを data.frame または matrix に既に読み込んでいるかどうかが明確ではありません。したがって、次の例では、テキスト ファイルにデータ セットがあると仮定します。

# read in data
dataset = read.table("dataset.txt", header=T)

# ... or use example data
dataset = read.table(text="ID X Y
1 118 10
2 119 50
3 100 12
4 101 12", header=T, sep=" ")

# create a matrix from columns X & Y and use project as in the question
project(as.matrix(dataset[,c("X","Y")]), "+proj=utm +zone=51 ellps=WGS84")
#             [,1]    [,2]
# [1,]   -48636.65 1109577
# [2,]   213372.05 5546301
# ...

アップデート：

コメントは、問題がproject()data.frame への適用に起因することを示唆しています。project()をチェックするため、data.frames では機能しませんis.numeric()。したがって、上記の例のようにデータを行列に変換する必要があります。を使用するコードに固執したい場合cbind()は、次のことを行う必要があります。

 X <- dd[,"X"]
 Y <- dd[,"Y"]
 xy <- cbind(X,Y)

dd["X"]との違いdd[,"X"]は、後者は data.frame を返さず、その結果cbind()、data.frame ではなくマトリックスを生成することです。

于 2013-09-05T15:40:20.647 に答える