2013-08-18 43 views
0

我從1983年開始在國家平面座標系統中具有以下數據集。所有座標位於長島區(3104)。將座標從長島區轉換爲經緯度。長格式與R

dput(example) 
structure(c(1008031L, 1000852L, 1001869L, 1005306L, 986887L, 
998031L, 1018703L, 1014319L, 1016186L, 1006977L, 1006891L, 1000883L, 
1001403L, 999812L, 1010077L, 1015918L, 984241L, 1013735L, 986848L, 
998243L, 1007312L, 1005663L, 992415L, 999771L, 1006787L, 987215L, 
990271L, 1015773L, 999342L, 1007245L, 1007098L, 996980L, 1006886L, 
999643L, 1008769L, 1016489L, 1004212L, 986848L, 1001512L, 1002584L, 
1001753L, 1004625L, 990725L, 1013435L, 1010795L, 1007509L, 1009419L, 
NA, 1009731L, 999007L, 999007L, 1000195L, 985863L, 990064L, 1008192L, 
1008306L, NA, 1003280L, 1006541L, 1001264L, 1003844L, 1008345L, 
987951L, 999104L, 1009013L, 998201L, 984182L, 1004940L, 1004513L, 
999659L, 1018204L, 1005918L, 1008158L, 999629L, 982208L, 1008008L, 
983985L, 1003591L, 992033L, 1012144L, 1008285L, 1004196L, 999937L, 
1007579L, 1001610L, 1013897L, 985504L, 1003588L, 1000088L, 1002230L, 
999304L, 1001393L, 997666L, 999148L, 997501L, 1004670L, 994699L, 
1005950L, 994821L, 998160L, 233036L, 228179L, 190702L, 186668L, 
173599L, 234924L, 241414L, 182198L, 178657L, 178140L, 242280L, 
236356L, 235184L, 238138L, 181374L, 245648L, 149582L, 211309L, 
212883L, 176387L, 243183L, 237170L, 149315L, 188471L, 242047L, 
215403L, 203844L, 240835L, 233575L, 234932L, 166665L, 174885L, 
193881L, 228852L, 244547L, 247336L, 178750L, 212883L, 232231L, 
248715L, 182080L, 242885L, 204176L, 251857L, 183147L, 245160L, 
235573L, NA, 243613L, 229814L, 229814L, 229856L, 212233L, 225331L, 
245037L, 245316L, NA, 229886L, 243541L, 232832L, 250988L, 235949L, 
220453L, 192913L, 242619L, 173610L, 150037L, 169914L, 180635L, 
229932L, 239783L, 190990L, 244973L, 243379L, 170319L, 246638L, 
205857L, 242274L, 215119L, 236944L, 245256L, 183865L, 238365L, 
183413L, 241367L, 238753L, 216029L, 249617L, 230093L, 176647L, 
227192L, 200533L, 177016L, 187285L, 170971L, 233870L, 176744L, 
179753L, 177866L, 227234L), .Dim = c(100L, 2L), .Dimnames = list(
    NULL, c("xcoord", "ycoord"))) 

我需要在lat./long有座標。格式或wgs84。

有人能告訴我我是怎麼做到的? 感謝您的幫助。

+0

這看起來有所幫助:http://gis.stackexchange.com/questions/16/how-do-i-convert-state-plane-coordinates-to-latitude-經度 – Thomas

+0

和這個http://cran.r-project.org/web/packages/proj4/proj4.pdf – yosukesabai

+0

感謝您的鏈接。我已經看到了兩個,但無法理解他們...... – Dominik

回答

3

替代rgdal的是proj4包。它具有transform的功能。我猜正在使用的投影是http://www.spatialreference.org/ref/epsg/32118/

library(proj4) 
dat <- structure(c(1008031L, 1000852L, 1001869L, 1005306L, 986887L, 
998031L, 1018703L, 1014319L, 1016186L, 1006977L, 1006891L, 1000883L, 
1001403L, 999812L, 1010077L, 1015918L, 984241L, 1013735L, 986848L, 
998243L, 1007312L, 1005663L, 992415L, 999771L, 1006787L, 987215L, 
990271L, 1015773L, 999342L, 1007245L, 1007098L, 996980L, 1006886L, 
999643L, 1008769L, 1016489L, 1004212L, 986848L, 1001512L, 1002584L, 
1001753L, 1004625L, 990725L, 1013435L, 1010795L, 1007509L, 1009419L, 
NA, 1009731L, 999007L, 999007L, 1000195L, 985863L, 990064L, 1008192L, 
1008306L, NA, 1003280L, 1006541L, 1001264L, 1003844L, 1008345L, 
987951L, 999104L, 1009013L, 998201L, 984182L, 1004940L, 1004513L, 
999659L, 1018204L, 1005918L, 1008158L, 999629L, 982208L, 1008008L, 
983985L, 1003591L, 992033L, 1012144L, 1008285L, 1004196L, 999937L, 
1007579L, 1001610L, 1013897L, 985504L, 1003588L, 1000088L, 1002230L, 
999304L, 1001393L, 997666L, 999148L, 997501L, 1004670L, 994699L, 
1005950L, 994821L, 998160L, 233036L, 228179L, 190702L, 186668L, 
173599L, 234924L, 241414L, 182198L, 178657L, 178140L, 242280L, 
236356L, 235184L, 238138L, 181374L, 245648L, 149582L, 211309L, 
212883L, 176387L, 243183L, 237170L, 149315L, 188471L, 242047L, 
215403L, 203844L, 240835L, 233575L, 234932L, 166665L, 174885L, 
193881L, 228852L, 244547L, 247336L, 178750L, 212883L, 232231L, 
248715L, 182080L, 242885L, 204176L, 251857L, 183147L, 245160L, 
235573L, NA, 243613L, 229814L, 229814L, 229856L, 212233L, 225331L, 
245037L, 245316L, NA, 229886L, 243541L, 232832L, 250988L, 235949L, 
220453L, 192913L, 242619L, 173610L, 150037L, 169914L, 180635L, 
229932L, 239783L, 190990L, 244973L, 243379L, 170319L, 246638L, 
205857L, 242274L, 215119L, 236944L, 245256L, 183865L, 238365L, 
183413L, 241367L, 238753L, 216029L, 249617L, 230093L, 176647L, 
227192L, 200533L, 177016L, 187285L, 170971L, 233870L, 176744L, 
179753L, 177866L, 227234L), .Dim = c(100L, 2L), .Dimnames = list(
    NULL, c("xcoord", "ycoord"))) 

lonlat <- project(dat, 
    '+proj=lcc +lat_1=41.03333333333333 +lat_2=40.66666666666666 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs ', 
    inverse=T) 
+0

不錯。最後,我認爲這兩種解決方案都使用PROJ4庫來進行投影。 – nograpes

+0

感謝您的回答。我通過Google搜索猜測了相同的EPSG投影。但你怎麼知道它是哪一個?你怎麼知道在項目命令中的''之間寫什麼?這樣我就不必再問下次是否有不同的投影...... ;-) – Dominik

+0

@Dominik,那些最好發佈到gis.stackexchange.com。最終,只有那些進行投影的人肯定知道,往往他們不知道確切的時間......我經常使用的一種方法是找到幾個可以在谷歌地球上識別的對象,並查看我的lon/lat反投影顯示出來。並確定錯誤的大小是我的應用程序可以接受的。 – yosukesabai

3

您可以使用rgdalsp來完成此操作。

library(rgdal) 
library(sp) 
# I made a guess at the PROJ4 string that your data was in. This is using SRID: 2831 
long.island.proj4<-CRS("+proj=lcc +lat_1=41.03333333333333 +lat_2=40.66666666666666 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs") 
# Remove the missing values in your data, and convert to SpatialPoints 
sp.points<-SpatialPoints(example[complete.cases(example),],proj4string=long.island.proj4) 
# Project to lat/long 
spTransform(sp.points,CRS('+proj=longlat')) 
+0

你是怎麼知道這是SRID:2831的?我在哪裏可以找到信息,CRS命令需要看起來像什麼?這樣我就不用再次詢問是否有不同的投影...... ;-) – Dominik

+1

@Dominik使用[在其他答案中提供的鏈接](http://www.spatialreference.org/ref/ EPSG /)。搜索「長島」,然後單擊PROJ4字符串。 – nograpes

+0

感謝您的回答 – Dominik