1) 經(jīng)緯度轉(zhuǎn)換
# DMS to DD
# method 1
> library(sp)
Warning message:
程輯包‘sp’是用R版本4.0.5 來建造的
> as.numeric(char2dms("48°26'5\"N", chd = "°", chm = "'", chs='"'))
[1] 48.43472
# method 2
# load test data
> head(test)
site lat lon
1 1 18°21'5" 43°59'5"
2 2 42°26'6" 113°30'6"
3 3 16°26'7" 97°34'7"
4 4 55°26'8" 99°13'8"
5 5 13°26'9" 20°42'9"
6 6 58°26'10" 109°23'10"
# extracting values from degree, minute, second with regular expression
> test1 = test %>%
mutate(Dlat = str_extract(test$lat, "(?<=^)(.*)(?=°)"),
Mlat = str_extract(test$lat, "(?<=°)(.*)(?=')"),
Slat = str_extract(test$lat, "(?<=')(.*)(?=\")"),
Dlon = str_extract(test$lon, "(?<=^)(.*)(?=°)"),
Mlon = str_extract(test$lon, "(?<=°)(.*)(?=')"),
Slon = str_extract(test$lon, "(?<=')(.*)(?=\")")) %>%
select(site, Dlat, Mlat, Slat, Dlon, Mlon, Slon) %>%
apply(., 2, as.numeric) %>%
as.data.frame() %>%
mutate(DDlat = Dlat + Mlat/60 + Slat/3600,
DDlon = Dlon + Mlon/60 + Slon/3600)
> head(test1)
site Dlat Mlat Slat Dlon Mlon Slon DDlat DDlon
1 1 18 21 5 43 59 5 18.35139 43.98472
2 2 42 26 6 113 30 6 42.43500 113.50167
3 3 16 26 7 97 34 7 16.43528 97.56861
4 4 55 26 8 99 13 8 55.43556 99.21889
5 5 13 26 9 20 42 9 13.43583 20.70250
6 6 58 26 10 109 23 10 58.43611 109.38611
2）樣點分布圖
ArcGIS也可以做，依然參照圖層概念著拭，在底圖上添加樣點纱扭，設置經(jīng)緯度網(wǎng)格等等。R中做法如下：
# load .shp files
# the map of china
> china <- read_sf("E:/GIS/四百萬數(shù)據(jù)/省界（線）.shp")
# method 1：tmap包
> library(sf)
> library(tmap)
> library(tmaptools)
> library(sp)
> library(readxl)
> library(grid)
# load information of site locations
> head(sites)
newname longitude latitude
1 Leo_1_R1 95.25550 31.54113
2 Leo_1_R2 95.25550 31.54113
3 Leo_1_R3 95.25550 31.54113
4 Leo_2_R1 92.93162 31.84779
5 Leo_2_R2 92.93162 31.84779
6 Leo_2_R3 92.93162 31.84779
# set spatial coordinates to create a Spatial object
coordinates(sites) <- ~longitude+latitude
proj4string(sites) <- CRS("+proj=longlat +datum=WGS84")
# plotting
> map1 <-
tm_shape(china,ylim=c(25,36), xlim=c(75,100))+
tm_lines(col = "darkgrey", lwd = 1)+
tm_layout(scale = 0.6,
legend.position = c("left","bottom"),
inner.margins = c(0.12,0.03,0.08,0.03))+
tm_shape(sites)+
tm_symbols(
size = 0.7,
shape = 16,
col = "royalblue3")+
tm_graticules(x= c(75,85,95,100),
y= c(25,28,31,34,37),
lwd = 0.5, col="lightgrey",
labels.size = 1.2)+
tm_compass(type = "arrow",
position = c("left","top"),
size = 5)+
tm_scale_bar(position = c("left","bottom"),
size = 1.2); map1
# method2: geom_sf()
# load information of climate stations in Tibet
> head(MAPsite)
# A tibble: 6 x 2
longitude latitude
<dbl> <dbl>
1 80.0 32.3
2 84.2 32.1
3 92.0 31.3
4 81.2 30.2
5 91.1 30.3
6 87.4 29.0
# convert foreign object to an sf object
> MAPsite = MAPsite %>%
st_as_sf(coords = c("longitude", "latitude")) %>%
st_set_crs(st_crs(china))
# plotting
> ggplot() +
geom_sf(data = china)+
geom_sf(data = MAPsite, col = "red")+
theme_bw()
Method 1
Method 2
