16EAGLE / basemaps

A lightweight package for accessing basemaps from open sources in R 🗺️
https://jakob.schwalb-willmann.de/basemaps
GNU General Public License v3.0
57 stars 15 forks source link

Add object to a basemap #6

Closed elgabbas closed 2 years ago

elgabbas commented 2 years ago

Hello,

I am interested in adding a spatial object to a base map (world_terrain_base) downloaded via the basemaps R package. I assume that the projection of the plotted base map object is epsg:3857, then I projected the object I want to plot into the same projection but was not able to overlay both together.

require(sf)
require(dplyr)
require(basemaps)
require(cleangeo)

set_defaults(map_service = "esri", map_type = "world_terrain_base")

# Area of interest (sf)
Area_sf <- rworldmap::getMap() %>% 
  cleangeo::clgeo_Clean() %>% 
  st_as_sf() %>% filter(ADMIN == "Egypt")
# Area of interest (sf): epsg:3857 projection
Area_sf_Prj <- st_transform(Area, 3857)

# Area of interest (sp)
Area_sp <- as_Spatial(Area_sf)
# Area of interest (sp): epsg:3857 projection
Area_sp_Prj <- spTransform(Area_sp, CRS("+init=epsg:3857"))

basemap_plot(Area_sf_Prj, add = T) # this works
plot(Area_sf_Prj[1], add = T) # this is not plotted
plot(Area_sp_Prj, add = T) # this is also not plotted

I also tried using the locator(1) function to know the range of coordinates of the plotted base map but it gives an unexpected range of coordinates.

Could you please help me to plot custom maps of interest to the base map? Is there another function from the package I should use instead of basemap_plot?

Thanks, Ahmed

asoultan commented 2 years ago

Dear Ahmed, The trick is with the projection. Basemap transforms the projection to another one. So to overlay a polygon or any other spatial formate, you need to to transform its projection first as follows:

require(sf) require(dplyr) require(basemaps) library(ggplot2) set_defaults(map_service = "esri", map_type = "world_terrain_base")

Area of interest as sf object

BM <- rworldmap::getMap() %>% cleangeo::clgeo_Clean() %>% st_as_sf() %>% filter(ADMIN %in% c("Egypt","Saudi Arabia"))# select two countries

Overlay layer

Egy <- rworldmap::getMap() %>% cleangeo::clgeo_Clean() %>% st_as_sf() %>% filter(ADMIN == c("Egypt")) Egypt<-as(Egy,"Spatial")

r<-basemap_raster(BM)# use this function to get the raster layer with 3 bands raster::projection(r)# to get to know the new projection

Egypt<-spTransform(Egypt,CRS("+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"))# tranform the overlay layer's projection to the new one to be consistant with basemap projection Neighborhoods <- fortify(Egypt)# to be overlaied by ggplot

s<-basemap_ggplot(BM)+ theme(axis.text.y = element_blank(), axis.text.x = element_blank(), axis.ticks = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank())+ xlab("") + ylab("")+ geom_polygon(aes(x=long, y=lat, group=group), fill='orange', size=.1,color='black', data=Neighborhoods, alpha=0.4) s# the final ggplot object

elgabbas commented 2 years ago

Thanks @asoultan. This works for me.

Cheers, Ahmed

16EAGLE commented 2 years ago

See also #7