[R-sig-Geo] Shapefile and Basemap

[sownal chand]

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.

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
>
>

	[[alternative HTML version deleted]]