r语言 - 如何将UTM转换为多头/拉特,同时点分布在不同的区域38,39,40,41上



我将把一万个UTM数据点转换为长/纬度协调系统。我的积分分布在伊朗的四个区域38,39,40和41。为了重现我的代码,我在以下代码中随机制作了 100 个样本集。如以下代码所示,当我选择区域 39 时,所有点都将显示在 leflet 地图上的区域 39 中。将区域更改为 38 时,所有点都将移动到区域 38,依此类推。我的问题是我如何处理此代码以将每个点协调到其真正正确的地理点。有没有办法从UTM中查找区域编号以使代码对不同区域有条件。谢谢

utmx<-c(453199, 452226, 459070, 456923 ,406917, 439882, 406917, 453199, 452226, 459070, 456923 ,439882, 452226 , 459070, 456923 ,406917, 439882 ,453199 ,452226 ,459070 ,456923 ,406917, 360103, 360103, 439000, 453400,414114, 413937, 410250, 402673 ,348019 ,346225 ,351019, 339701, 335818, 361243, 360250, 379132, 375932 , 379140, 362556, 360103, 423200, 439000, 414114, 413937, 410250, 402673, 393698, 361243, 356800, 365150 , 360250, 379132, 375932, 379140, 393917, 401315, 362556, 360103, 348019, 346225, 351019, 335818, 338793, 393698, 361243, 356800, 365150, 360250, 362556, 360103, 439000, 414114, 413937, 410250, 402673, 348019 , 351019, 339701, 335818, 333066, 338793, 326138, 361243, 356800, 365150, 360250, 379132, 379140, 393917 , 401315, 360103, 414114, 413937, 410250, 402673, 393698, 356800, 365150)

utmy<-c(4094292, 4095754, 4091838, 4092360, 4084991, 4095385 ,4084991, 4094292, 4095754, 4091838, 4092360, 4095385, 4095754, 4091838, 4092360 ,4084991, 4095385 ,4094292, 4095754 ,4091838 ,4092360, 4084991,4067970 ,4067970 ,4020600 ,4020400, 4081856 ,4083323 ,4082400, 4087420, 4102943 ,4115011, 4111792, 4139663, 4133857, 4071241, 4078450, 4059477, 4059055, 4059506 ,4069070 ,4067970, 4113100, 4020600, 4081856, 4083323, 4082400 ,4087420, 4068252, 4071241, 4082850, 4078100, 4078450,4059477 ,4059055, 4059506, 4063361, 4069043, 4069070, 4067970, 4102943, 4115011 ,4111792, 4133857 ,4120862, 4068252,4071241, 4082850 ,4078100 ,4078450, 4069070 ,4067970 ,4020600 ,4081856 ,4083323 ,4082400, 4087420,4102943, 4111792, 4139663, 4133857, 4138934 ,4120862 ,4131102 ,4071241, 4082850, 4078100, 4078450, 4059477, 4059506, 4063361, 4069043, 4067970, 4081856, 4083323 ,4082400 ,4087420 ,4068252,4082850,
4078100)
body<-data.frame(utmx,utmy)
###For Iran
wgs84 = "+init=epsg:4326"  #this need to be checked
bng = "+proj=utm +zone=39 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
### 
##########Conversion function
library(rgdal)
ConvertCoordinates <- function(easting,northing) {
out = cbind(easting,northing)
mask = !is.na(easting)
sp <-  sp::spTransform(sp::SpatialPoints(list(easting[mask],northing[mask]),proj4string=sp::CRS(bng)),sp::CRS(wgs84))
out[mask,]=sp@coords
out
}
########### Conversion using above function and checking with web convertor
LongLat<-ConvertCoordinates(body$utmx,body$utmy)
colnames(LongLat)<-c("lon","lat")
body<-cbind(body, LongLat ) 
###############################first prepare a leaflet plot ...
lplot <- leaflet::leaflet(data = body) %>% # create leaflet object
leaflet::addTiles() %>% # add basemap
#leaflet::addCircleMarkers( clusterOptions = markerClusterOptions())# add data layer - markers
leaflet::addMarkers(clusterOptions = markerClusterOptions())
lplot #  ... then display it

如果您只有坐标号而不是区域标识符,则无法知道您的点位于哪个UTM区域中。UTM 坐标号对每个区域重复,只有区域标识符才能告诉您点所在的区域。点(453199,4094292(可能在世界任何地方,更不用说伊朗的这四个区域了。

您需要返回到源并获取每个点的UTM区域信息,然后您可以为每个点构造正确的投影字符串,并为每个不同的UTM区域进行变换,以获得经度坐标。

最新更新