Shapefile and Basemap

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 at 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 at coords)
    >
    > head(data.df_SPDF at 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