[R-sig-Geo] Shapefile and Basemap

[Micha Silver]

On 5/15/22 03:52, sownal chand wrote:
> Hi Micha,
>
> Thanks for your email and really appreciate your assistance.
>
> I see that there are alot of old versions of information on the net 
> with spatial analysis using R. I would like to request if you have any 
> particular site or reading materials for spatial analysis and 
> modelling in regards to R programming would really help.
>

One of the best online sources (in my opinion) is the book by Robin 
Lovelace "Geocomputation with R":

https://geocompr.robinlovelace.net/


Here's a good explanation of all the changes in the R-spatial packages, 
and why those changes were necessary:

https://r-spatial.org/r/2022/04/12/evolution.html


If you're following tutorials that depend on any of the packages: `sp`, 
`raster`, `rgdal`, then you should be aware that these will soon be 
deprecated and your procedure will not work in a year or so.

(Authors of those tutorials will have to revise them...)


Good luck



> Once again, thanking you in advance
> Kind regards
> Sownalc
>
> On Sun, May 15, 2022, 03:50 Micha Silver <tsvibar using gmail.com> wrote:
>
>     Hello
>
>
>     On 5/13/22 12:47, sownal chand wrote:
>     > Hello sir/madam,
>     >
>     > I am working with shape file of my country and the issue I am
>     facing
>     > is the shapefile is scattered while plotting it using basemap. I am
>     > using my sample point data which is attached to this email. I hope
>     > that some expert in this area would help in correcting the codes
>     below
>     > to show the shapefile in one location ( its a pacific centered map)
>
>
>     When you raised this question a few weeks ago, it was suggested to
>     avoid
>     the `sp` package with its SPDF data type, and instead focus on the
>     newer
>     `sf` package.
>
>     There is also a replacement for the raster::getData function in
>     the new
>     `geodata` package.
>
>
>     Here is a much simpler version of what (I think) you are trying to
>     achieve:
>
>
>     # Load only three libraries
>
>     library(sf)
>     library(tmap)
>     library(geodata)
>
>
>     # Read your list of data (You should remove the summary line in
>     advance...)
>
>     dt <- read.csv("DataR.csv")
>     dt <- dt[complete.cases(dt),]
>     dt_sf <- st_as_sf(dt, coords=c("long", "lat"),
>                        crs="EPSG:4326")
>     str(dt_sf)
>
>     # Get Fiji boundary from geodata package
>     fiji <- gadm(country="FIJI", level=2, path=tempdir())
>     # Convert to sf object for tmap plotting
>     fiji <- st_as_sf(fiji)
>
>
>     # Visualize with tmap
>
>     tmap_mode("view")
>     tm_basemap("OpenStreetMap.Mapnik") +
>              tm_shape(fiji) +
>              tm_borders(col="brown", lwd=2) +
>              tm_shape(dt_sf) +
>
>              # size of symbols by yearly data. You can choose any
>     year, of
>     course
>
>              tm_symbols(col="blue", size="Year2", scale=0.1)
>
>
>
>     If you have a specific problem with a certain shapefile, you'll
>     have to
>     supply it to the list in order to get further help.
>
>
>     HTH
>
>
>     >
>     >
>     **************************************************************************************************
>     >
>     >
>     > library(sp)
>     > library(raster)
>     > library(rgdal)
>     > library(leaflet)
>     >
>     > #read.csv
>     > read.csv ("C://Users/Documents/data.csv") -> data.df
>     > head(data.df)
>     >
>     > hist(data.df$Year, breaks=20)
>     >
>     > #remove NA valuues in the spatial Data Frame
>     > data.df <- na.omit(data.df)
>     > View(data.df)
>     >
>     > plot(data.df$long, data.df$lat,
>     >      ylab = "Latitude", xlab="Longitude") #boring!
>     >
>     > # Use the cex function to plot circle size as a function of a
>     variable
>     > plot(data.df$long, data.df$lat,
>     >      cex = data.df$Year.7 * 0.045,
>     >      ylab = "Latitude", xlab="Longitude")
>     >
>     > data.df_SPDF <- SpatialPointsDataFrame(coords = data.df[,c("long",
>     > "lat")],
>     >                                                 data =
>     > data.df[,c("Year", "Year.1", "Year.2","Year.3","Year.4")],
>     > proj4string =
>     > CRS("+init=epsg:4326")) # sets the projection to WGS 1984 using
>     > lat/long. Optional but good to specify
>     >
>     > # Summary of object
>     > data.df_SPDF
>     >
>     > # SPDFs partition data elements, e.g. the coordinates are stored
>     > separately from the data
>     > head(data.df_SPDF using coords)
>     >
>     > head(data.df_SPDF using data)
>     >
>     >
>     > # You can quickly access the data frame as per a standard data
>     frame, e.g.
>     > head(data.df_SPDF$Year)
>     >
>     >
>     > # You can use the plot or spplot function to get quick plots
>     > plot(data.df_SPDF)
>     >
>     > spplot(data.df_SPDF, zcol = "Year")
>     >
>     > FIJ_Adm_2 <- readOGR("FIJ_Adm_2_shapefile", "FIJ_Adm_2")
>     >
>     > # You first need the ISO3 codes for the country of interest.
>     >
>     > # The getData function then allows you to retrieve the relevant
>     admin
>     > level boundaries from GADM.
>     > FJI_Adm_2 <- raster::getData("GADM", country="FJI", level=2)
>     >
>     > # Plot both country and data points
>     > plot(FJI_Adm_2)
>     > points(data.df$long, data.df$lat,
>     >        cex = data.df$Year * 0.045,
>     >        ylab = "Latitude", xlab="Longitude",
>     >        col="red")
>     >
>     > basemap <- leaflet() %>% addProviderTiles("CartoDB.Positron")
>     >
>     > basemap %>% addPolygons(data=FJI_Adm_2)
>     >
>     > # to change the colors/line weight
>     > basemap %>% addPolygons(data=FJI_Adm_2, color = "red",
>     >                         weight = 1, fillOpacity = 0.2)
>     >
>     >
>     > # If you want to add points as well
>     > basemap %>% addPolygons(data=FJI_Adm_2, weight = 2,
>     >                         popup = FJI_Adm_2$NAME_2) %>%
>     >
>     >   addCircleMarkers(data=data.df_SPDF,
>     >                    color="red", radius = 2)
>     >
>     > library(wesanderson) # for a nice color palette
>     > colorPal <- colorNumeric(wes_palette("Zissou1")[1:5],
>     > data.df_SPDF$Year, n = 5)
>     >
>     >
>     > # colorPal is now a function you can apply to get the corresponding
>     > # color for a value
>     > colorPal(0.6)
>     >
>     > basemap %>% addPolygons(data=FJI_Adm_1, weight = 2, fillOpacity=0,
>     >                         popup = FJI_Adm_1$NAME_1) %>%
>     >
>     >   addCircleMarkers(data=data.df_SPDF,
>     >                    color = colorPal(data.df$Year),
>     >                    radius = 2,
>     >                    popup = as.character(data.df$Year))%>%
>     >
>     >   addLegend(pal = colorPal,
>     >             title = "Average Temp for Year",
>     >             values = data.df_SPDF$Year)
>     >
>     >
>     >
>     >
>     *************************************************************************************************
>     >
>     > Thanking you in advance
>     > sownalc
>
>     -- 
>     Micha Silver
>     Ben Gurion Univ.
>     Sde Boker, Remote Sensing Lab
>     cell: +972-523-665918
>     https://orcid.org/0000-0002-1128-1325
>
-- 
Micha Silver
Ben Gurion Univ.
Sde Boker, Remote Sensing Lab
cell: +972-523-665918
https://orcid.org/0000-0002-1128-1325