我无法找到多边形和空间点之间的距离。我试图测量区域(多边形)和以前的锡矿(空间点)之间的距离。我怀疑我可能在投影方面做错了。我似乎能够计算出各区之间的距离,但不能计算区和矿的距离。
#Open file
mapping<- readOGR(dsn="C:/Users/noble/OneDrive - London School of Economics/LSE/EC465/Extended Essay/Stack Exchange Question/gadm2.Peninsula.dbf")
#Define the districts
Kinta <- mapping[mapping@data$NAME_1 == "Perak" &
mapping@data$NAME_2 == "Kinta", ]
Gombak <- mapping[mapping@data$NAME_1 == "Selangor" &
mapping@data$NAME_2 == "Gombak", ]
#Set a projection
EPSG.3375<-"+proj=omerc +lat_0=4 +lonc=102.25 +alpha=323.0257964666666 +k=0.99984 +x_0=804671 +y_0=0 +ellps=GRS80 +units=m +no_defs"
Gombak.km<-spTransform(Gombak,CRS(EPSG.3375))
Kinta.km<-spTransform(Kinta,CRS(EPSG.3375))
#Calculate the distance
gDistance(Kinta.km, Gombak.km, byid=TRUE)
#Open data on tin mines and define as spatial points
tin.01<-read.csv("C:/Users/noble/OneDrive - London School of Economics/LSE/EC465/Extended Essay/Stack Exchange Question/Tin Mine Location_01.csv")
coordinates(tin.01)<-8:9
proj4string(tin.01) <- CRS("+proj=omerc +lat_0=4 +lonc=102.25 +alpha=323.0257964666666 +k=0.99984 +x_0=804671 +y_0=0 +ellps=GRS80 +units=m +no_defs")
#Find distance between district and mines
gDistance(Kinta.km,tin.01,byid=TRUE)
多边形之间距离的输出似乎是正确的:
> gDistance(Kinta.km, Gombak.km, byid=TRUE)
59
71 100676
但我对地区和矿山之间距离的输出肯定是错误的:
> gDistance(Kinta.km,tin.01,byid=TRUE)
59
1 661153.5
2 661152.6
3 661153.0
4 661152.7
5 661151.8
6 661152.9
7 661153.1
8 661153.3
9 661153.2
我打算做的是:1)计算马来西亚所有地区与所有以前的锡矿之间的距离; 2)提取最近的锡矿到每个区的距离。
这是我正在使用的数据的link。我正在使用区域多边形的GADM数据,并且我亲自编写了历史锡矿的位置。非常感谢所有帮助。谢谢。
您当前的方法存在两个问题。 1.)coordinates
任务不太正确。它应该是long / lat,但是你将它指定为lat / long。 2.)以当前方式直接设置CRS实际上不会以必要的方式改变点。您需要首先分配适当的长/纬度CRS,然后执行spTransform
操作。
#Open data on tin mines and define as spatial points
tin.01 <- read.csv("random/Tin Mine Location_01.csv")
coordinates(tin.01) <- c("Longitude","Latitude") ## Should be `9:8` if you wanted to stick with indexes, but using the names here is generally lower risk.
proj4string(tin.01) <- CRS(proj4string(Kinta)) ## Setting the initial projection to the same one the polygons are using. You should change this if your original data source uses some other known long/lat projection.
tin.km <- spTransform(tin.01,CRS(EPSG.3375)) ## Creating a transformed set of points for the distance calculation.
#Find distance between district and mines
gDistance(Kinta.km,tin.km,byid=TRUE) ## '0' distance means the mine is inside the district.
59
1 194384.372
2 223773.999
3 0.000
4 36649.914
5 102944.361
6 0.000
7 0.000
8 6246.066
9 0.000