4

我在这里绘制地图，无法正确填写，希望您能提供帮助。那我该怎么办

阅读地图并更改投影：

map<-readOGR(dsn="e:\\r\\OSM_adm_reg\\adm4_region (1)", layer="adm4_region")
map=spTransform(map,CRS("+proj=latlong +ellips=GRS80"))

读取全国各地物流仓库位置数据，匹配world.cities数据库在地图上绘制：

    data(world.cities)
    russian.cities=world.cities[world.cities$country.etc=="Russia",]
    metro_outlets=read.csv(file="e:\\r\\retail\\metro.csv",head=TRUE,sep=";")
    match<-match(metro_outlets$CITY_ENGLISH, russian.cities$name)
    outlet_coordinates=russian.cities[match,]
    outlet_coordinates=cbind(outlet_coordinates, "metro_store_count"=metro_outlets$STORE_COUNT_METRO)
    outlet_coordinates=cbind(outlet_coordinates, "cch_depo"=metro_outlets$DEPO_FOOD)
    outlet_coordinates=cbind(outlet_coordinates, "NAME_LAT"=metro_outlets$STORE_OBLAST)
    outlet_coordinates=cbind(outlet_coordinates, "cch_franchise"=metro_outlets$Franchise)
outlet_coordinates=cbind(outlet_coordinates, "SL"=metro_outlets$SL)

基本上，我需要根据 SL 列的值填充区域。现在，我试图获取包含 WH 的区域列表：

outlet_spatial=SpatialPoints(outlet_coordinates[,c(5,4)],proj4string=CRS(proj4string(map)))
index=over(outlet_spatial,map)
region_names=rownames(map@data[map@data$NAME_LAT %in% index$NAME_LAT,])
map_df<-fortify(map,map$NAME_LAT)
map_df$fill=ifelse(map_df$id %in% region_names,"true","false")

并绘制地图：

ggplot()+
  geom_polygon(data=map_df,aes(x=long,y=lat,group=group,fill=fill),col="gray")+
  scale_fill_manual(values=c("gray","gray70"))+xlim(20, 180) + ylim(40,80)+
  geom_point(aes(x=long,y=lat), color="orange",alpha=I(8/10), data=outlet_coordinates)

这是我得到的地图！

现在，我不知道如何将 SL 列的值传递给map_df并相应地填充区域。SL 在这里定义： outlet_coordinates=cbind(outlet_coordinates,"SL"=metro_outlets$SL)

4

1 回答 1

3

这比我想象的要复杂一些，主要是因为你得到的 TOPOLOGY 错误是由你的 shapefile 中糟糕的几何形状引起的。也许您可以改用高质量、公开可用的 shapefile。我建议www.naturaleartchdata.com。我在脚本中包含了一个链接并下载了该数据。您可以复制并粘贴此内容。这应该有效：

require(ggplot2)
require(maptools)
require(sp)

# Your shape file has some topology errors meaning it can't be used by 'over()'
# Use open source shapefiles from www.naturalearth.com - v. good quality
oldwd <- getwd(); tmp <- tempdir(); setwd(tmp)
url <- "http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_admin_1_states_provinces_shp.zip"
dest <- paste(tmp,"\\tmp.zip",sep="")
download.file( url , dest )
unzip(dest)


#   Projection information for geographic data (but in decimal degrees)
projgeo <- CRS("+proj=latlong +datum=WGS84 +ellps=WGS84")


# Read in world shapefile
wld <- readShapePoly("ne_10m_admin_1_states_provinces_shp" , proj4string = projgeo)
# Subset to Russia
map <- wld[wld$admin == "Russia" , ]


# City and metro data
r.city <- world.cities[ world.cities$country.etc=="Russia" , ]
metro_outlets <- read.csv(file="e:\\r\\retail\\metro.csv",head=TRUE,sep=";")


# Assign metro outlets to city and make spatial point data frame
met <- subset( metro_outlets , select = c( "SL" , "CITY_ENGLISH" ) )
met$SL <- as.numeric( met$SL )
match <- match( met$CITY_ENGLISH , r.city$name )
coordinates( met ) <- r.city[ match , c( "long" , "lat" ) ]
proj4string( met ) <- proj4string( map )


# Assign metro outlet attribute "SL" to each region of map
# Find the average SL for each region using 'over()'
df <- over( map , met[ , "SL" ] , fn = mean , na.rm = TRUE )
map$SL <- as.numeric( df$SL )


# Convert th map into format for plotting with ggplot
map@data$id <- rownames( map@data )
map.df <- as.data.frame( map )
map.fort <- fortify( map , region = "id" )
map.gg <- join( map.fort , map.df , by = "id" )


# Plot  
p1 <-   ggplot( NULL ) + 
    geom_polygon( mapping = aes( long , lat , group = group , fill  = factor( SL ) ) , colour = "white" , data = map.gg , map = map.gg )+
    scale_x_continuous( limits = c( 20 , 180 ) )
print(p1)

在此处输入图像描述

希望有帮助！

于 2013-03-09T23:33:18.810 回答