代码之家  ›  专栏  ›  技术社区  ›  sarovasta

映射点和多边形

  •  1
  • sarovasta  · 技术社区  · 9 年前

    我是一个R初学者,第一次处理R和空间数据。所以我希望我说清楚。

    我有一个意大利地区的形状文件,分为几个人口普查区。 另一方面,我有一个csv文件,其中包含案例列表和每个案例的地址。 我想映射点和形状文件,并获得每个人口普查分区中有多少个点的计数。

    #get cases file
    cases <- read.csv("cases.csv", sep =';', header = TRUE)
    names(cases)
    [1] "name"    "address"
    
    #geocode addresses from Google Maps
    library(GISTools)
    library(rgeos)
    library(ggmap)
    geolocalize <- geocode(as.character(cases$address))
    
    # bind latitude and longitude to the previous cases data frame
    cases <- data.frame(cases, geolocalize)
    names(cases)
    [1] "name"    "address" "lon"     "lat"
    
    #make cases a SpatialPointDataFrame
    #since addresses were retrieved using GoogleMaps, I set proj4string as follows
    cases.points <- SpatialPointsDataFrame(cases[,3:4], cases, proj4string = CRS("+init=EPSG:3857"))
    
    #get the shapefile
    region <- readOGR("R02_11_WGS84.shp")
    

    现在,我能够策划案件。点和形状文件分开,但不能将它们添加到同一绘图中。此外,正如我所说,我想计算每个多边形中有多少个点(即“区域”的人口普查划分)。

    我必须承认我对地理不是很感兴趣。我怀疑不同的坐标和/或投影参考系可能是问题所在,因此我检查了一下。

    head(coordinates(region))
    [,1]    [,2]
    0 364509.0 5065900
    1 363916.3 5056629
    2 372585.0 5068078
    3 360692.3 5048321
    4 356029.7 5062399
    5 360012.1 5065663
    
    
    coordinates(cases)
    lon      lat
    [1,] 7.323667 45.73664
    
    proj4string(region)
    [1] "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
    
    proj4string(cases.points)
    [1] "+init=EPSG:3857 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"
    

    形状文件坐标是否可能以度为单位,而案例坐标是否以十进制为单位?如果是,如何转换?

    萨罗

    1 回复  |  直到 9 年前
        1
  •  2
  •   user3757897 jdharrison    9 年前

    陈列 地图,但如果你下载纵横比,那只是未投影的坐标。要正确读取横向/纵向数据,请使用:

    library(sp)
    cases.points <- SpatialPointsDataFrame(cases[,3:4], cases, 
                                       proj4string = CRS("+init=EPSG:4326"))
    

    然后,对于投影,您需要spTransform-尝试:

    cases.points.utm32 <- spTransform(cases.points, CRS(proj4string(region)))
    

    编辑: 要在多边形内选择点,需要over()函数,也需要来自sp包(这是一个单独的问题)。一定要阅读sp包的基本函数和sp类的工作原理-下面基本上是从over()的帮助部分复制的。

    region$pointCount <- sapply(over(region, geometry(cases.points), 
                                 returnList = TRUE), length)
    
    推荐文章