以下是我拥有的数据框的示例,该数据框是从可公开获取的圣路易斯犯罪数据集中获得的。与数据相关的文档指出Xcoord和Ycoord位于[State Plane North America Datum 1983(NAD83)格式
CodedMonth Description XCoord YCoord
1: 2019-09 AUTO THEFT-PERM RETNT/UNRECOV OVER 48HR 908297.3 1018623.0
2: 2019-09 ASSLT-AGGRAV-OTH-WPN-2ND-CHILD-DOMESTIC 903995.7 1014255.0
3: 2019-09 FORGERY-ISSUING FALSE INSTRUMENT OR CERTIFICAT 0.0 0.0
4: 2019-09 STLG BY DECEIT/IDENTITY THEFT REPORT 890704.7 1010659.0
5: 2019-09 STALKING (HARASSMENT ONLY, NO THREAT) 881105.8 1008297.0
6: 2019-09 LARCENY-MTR VEH PARTS UNDER $500 882929.6 992941.3
如何将它们转换为Xcoord
和Ycoord
列为lon和lat格式,以便可以使用ggmap进行绘制
我找到了几个答案Convert latitude/longitude to state plane coordinates
但是我似乎无法使其用于我的数据
您可以使用sf库将其转换为简单的要素地理。为了使它起作用,您需要知道使用的坐标系,并根据提供的描述(国家平面NAD83,位于圣路易斯附近),我猜EPSG 26996(NAD83 /密苏里州东部USFT)是一个很好的起点。您可以在spatialreference.org处查找投影。更改为crs = ESRI:102696,因为第一个猜测无效。
library(sf)
# assuming your data is an object called my_df
my_sf_df <- st_as_sf(my_df, coords = c("XCoord", "YCoord"), crs = 102696) # 26996 didn't plot right
这会将x和y设置为空间坐标。您需要重新投影到WGS84这样的地理系统中,才能转换为纬度,您可以使用st_transform
my_latlon_df <- st_transform(my_sf_df, crs = 4326 ) %>%
mutate( lat= st_coordinates(my_latlon_df)[,1],
lon = st_coordinates(my_latlon_df)[,2])
my_latlon_df
# Simple feature collection with 6 features and 5 fields
# geometry type: POINT
# dimension: XY
# bbox: xmin: -93.26566 ymin: 35.80151 xmax: -90.19163 ymax: 38.63065
# epsg (SRID): 4326
# proj4string: +proj=longlat +datum=WGS84 +no_defs
# # A tibble: 6 x 6
# X1 CodedMonth Description geometry lat lon
# * <chr> <chr> <chr> <POINT [°]> <dbl> <dbl>
# 1 1: 2019-09 AUTO THEFT-PERM RETNT/UNRECOV OVER 48HR (-90.19163 38.63065) -82.2 44.7
# 2 2: 2019-09 ASSLT-AGGRAV-OTH-WPN-2ND-CHILD-DOMESTIC (-90.20674 38.6187) -82.3 44.7
# 3 3: 2019-09 FORGERY-ISSUING FALSE INSTRUMENT OR CERTIFICAT (-93.26566 35.80151) -93.3 35.8
# 4 4: 2019-09 STLG BY DECEIT/IDENTITY THEFT REPORT (-90.25329 38.60893) -82.4 44.6
# 5 5: 2019-09 STALKING (HARASSMENT ONLY, NO THREAT) (-90.2869 38.60251) -82.5 44.6
# 6 6: 2019-09 LARCENY-MTR VEH PARTS UNDER $500 (-90.28065 38.56034) -82.5 44.5
box <- st_bbox(my_latlon_df)
names(box) <- NULL
base_map <- get_map(location = box, zoom = 7)
ggmap(base_map)+
geom_sf(data = my_latlon_df,
aes(color = CodedMonth)
)
似乎没有坐标的点引起了问题,您需要确定如何处理。也想绕点]
# let's exclude point 3 for now
my_latlon_df <- my_latlon_df[c(1:2, 4:6),]
box <- st_bbox(my_latlon_df) # bounding box
names(box) <- NULL # removing non-complient labels
box2 <- box + c(-.2, -.2, .2, .2) # buffering
base_map <- get_map(location = box2, source = "osm") # getting base map
ggmap(base_map)+
geom_sf(data = my_latlon_df,
aes(color = CodedMonth)
)+
scale_x_continuous(limits = c(-90.3, -90.1))+
scale_y_continuous(limits = c(38.5, 38.7))