Reputation: 43
I'm plotting a map here and can't fill it right, I hope you can help. So what do I do
Read the map and change projection:
map<-readOGR(dsn="e:\\r\\OSM_adm_reg\\adm4_region (1)", layer="adm4_region")
map=spTransform(map,CRS("+proj=latlong +ellips=GRS80"))
Read data about location of logistics warehouses throughout the country and match the world.cities database to plot them on the map:
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)
Basically, I'm gonna need to fill regions based on the values of the SL column. Now, Im trying to get a list of regions that contains the WHs:
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")
And plot the map:
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)
This is the map I got!
Now, I don't know how to pass the values of the SL column to the map_df
and fill regions accordingly. SL defined here: outlet_coordinates=cbind(outlet_coordinates,"SL"=metro_outlets$SL)
Upvotes: 4
Views: 3518
Reputation: 59970
This was a little more involved than I thought, mainly because the TOPOLOGY errors you got are caused by poor geometry in your shapefile. Perhaps you could instead use good quality, publically available shapefiles. I suggest www.naturaleartchdata.com. I have included a link and downloaded that data in the script. You can copy and paste this. This should work:
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)
Hope that helps!
Upvotes: 3