CODEWITHSUNDEEP
rcoordinatesgeospatialspatialr-sf
BiggestNoob
BiggestNoob

Reputation: 1

How do I create a grid in R and plot my coordinates?

any help would be very appreciated. I am trying to plot some data I collected for a class project but am having trouble figuring it out, because, as the title suggests, I'm a complete noob.

I have 70 plot points of American Robin occurrence data throughout a burned area. My data set includes two separate visits with 1 or 0 values, (1 meant a Robin was observed, 0 was not observed), the plot ID for each point, and the X and Y coordinates for each of 70 points in Decimal Degrees. I want to plot these in R on a grid based on the shape of the wildfire raster I have.

I’ve tried converting my data frame to an sf object, but am very lost. Any help is appreciated. Previously, I've only had to work with data in the class that was already in the R environment, so it was much easier for me to get started. Thanks very much for any insight or direction.

I tried this, but am a beginner at R so unsurprisingly it does not work:

Robindf = as.data.frame(Robin)

# Convert data frame to sf object
RobinSpatialFrame <- st_as_sf(x = Robindf, 
                              coords = c("Lon", "Lat"))

Thank you!!

Here is what my Robindf looks like:

head(Robindf, 10) |> dput()

structure(list(PlotID = c("HS01", "HS02", "HS04", "HS05", "HS06", 
"HS07", "HS08", "HS09", "HS10", "HS11"), Lon = c(-105.1359899, 
-105.1481566, -105.297351, -105.2949066, -105.292101, -105.2905732, 
-105.2101288, -105.1492122, -105.1459622, -105.1448767), Lat = c(39.0979659, 
39.11979923, 39.23021589, 39.23077145, 39.23110478, 39.22882701, 
39.09652146, 39.12149368, 39.12157701, 51.08832563), `Visit 1` = c(0, 
1, 0, 0, 0, 0, 1, 1, 1, 1), `Visit 2` = c(0, 0, 0, 0, 0, 0, 0, 
0, 0, 0)), row.names = c(NA, 10L), class = "data.frame")

head(hayman, 10) |> dput()
 
 
structure(list(Event_ID = c("CO3922010528720020608", "CO3922010528720020608", 
"CO3922010528720020608", "CO3922010528720020608", "CO3922010528720020608", 
"CO3922010528720020608"), irwinID = c(NA_character_, NA_character_, 
NA_character_, NA_character_, NA_character_, NA_character_), 
    Incid_Name = c("HAYMAN", "HAYMAN", "HAYMAN", "HAYMAN", "HAYMAN", 
    "HAYMAN"), Incid_Type = c("Wildfire", "Wildfire", "Wildfire", 
    "Wildfire", "Wildfire", "Wildfire"), Map_ID = c("12585", 
    "12585", "12585", "12585", "12585", "12585"), Map_Prog = c("MTBS", 
    "MTBS", "MTBS", "MTBS", "MTBS", "MTBS"), Asmnt_Type = c("Extended", 
    "Extended", "Extended", "Extended", "Extended", "Extended"
    ), BurnBndAc = c("128726", "265", "322", "29", "50", "25"
    ), BurnBndLat = c("39.15", "39.316", "39.228", "39.224", 
    "39.153", "39.253"), BurnBndLon = c("-105.268", "-105.233", 
    "-105.254", "-105.171", "-105.438", "-105.355"), Ig_Date = c("2002/06/08", 
    "2002/06/08", "2002/06/08", "2002/06/08", "2002/06/08", "2002/06/08"
    ), Pre_ID = c("503303320010824", "503303320010824", "503303320010824", 
    "503303320010824", "503303320010824", "503303320010824"), 
    Post_ID = c("503303320030814", "503303320030814", "503303320030814", 
    "503303320030814", "503303320030814", "503303320030814"), 
    Perim_ID = c(NA_character_, NA_character_, NA_character_, 
    NA_character_, NA_character_, NA_character_), dNBR_offst = c("29", 
    "29", "29", "29", "29", "29"), dNBR_stdDv = c("-9999", "-9999", 
    "-9999", "-9999", "-9999", "-9999"), NoData_T = c("-970", 
    "-970", "-970", "-970", "-970", "-970"), IncGreen_T = c("-150", 
    "-150", "-150", "-150", "-150", "-150"), Low_T = c("140", 
    "140", "140", "140", "140", "140"), Mod_T = c("211", "211", 
    "211", "211", "211", "211"), High_T = c("350", "350", "350", 
    "350", "350", "350"), Comment = c(NA_character_, NA_character_, 
    NA_character_, NA_character_, NA_character_, NA_character_
    )), row.names = 0:5, class = "data.frame")

the CRS for the Hayman shapefile shows up as this:

crs(hayman)

Coordinate Reference System:
Deprecated Proj.4 representation:
 +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs 
WKT2 2019 representation:
PROJCRS["USA_Contiguous_Albers_Equal_Area_Conic_USGS_version",
    BASEGEOGCRS["NAD83",
        DATUM["North American Datum 1983",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["USA_Contiguous_Albers_Equal_Area_Conic_USGS_version",
        METHOD["Albers Equal Area",
            ID["EPSG",9822]],
        PARAMETER["Latitude of false origin",23,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",-96,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",29.5,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",45.5,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Not known."],
        AREA["United States (USA) - CONUS onshore - Alabama; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming."],
        BBOX[24.41,-124.79,49.38,-66.91]],
    ID["ESRI",102039]]

Upvotes: 0

Views: 830

Answers (1)

dimfalk
dimfalk

Reputation: 1413

In order to create a simple feature object, sf::st_as_sf() was already a good pick - just don't forget to specify your crs = "epsg:4326" argument to provide a spatial reference to your coordinates:

library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE

# import
p <- sf::st_as_sf(Robindf, coords = c("Lon", "Lat"), crs = "epsg:4326")
p
#> Simple feature collection with 10 features and 3 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -105.2974 ymin: 39.09652 xmax: -105.136 ymax: 51.08833
#> Geodetic CRS:  WGS 84
#>    PlotID Visit 1 Visit 2                   geometry
#> 1    HS01       0       0  POINT (-105.136 39.09797)
#> 2    HS02       1       0  POINT (-105.1482 39.1198)
#> 3    HS04       0       0 POINT (-105.2974 39.23022)
#> 4    HS05       0       0 POINT (-105.2949 39.23077)
#> 5    HS06       0       0  POINT (-105.2921 39.2311)
#> 6    HS07       0       0 POINT (-105.2906 39.22883)
#> 7    HS08       1       0 POINT (-105.2101 39.09652)
#> 8    HS09       1       0 POINT (-105.1492 39.12149)
#> 9    HS10       1       0  POINT (-105.146 39.12158)
#> 10   HS11       1       0 POINT (-105.1449 51.08833)

Let's assume this is correct if you're working in the area of Denver - feel free to plot your data via plot(p). If you wanted, you could also export your data as shapefile using sf::st_write(p, "robin.shp).

If you wanted to overlay your points with your hayman dataset in R, you have to specify the coordinate reference system of your vector dataset once again, and better also reproject your data to a common crs. The EPSG code 102039, as shown in your WKT2 notation, reveals ESRI being the issuing authority here.

# import and reproject to WGS 84, EPSG: 4326
grid <- sf::st_as_sf(hayman, coords = c("BurnBndLat", "BurnBndLon"), crs = "ESRI:102039") |>
  sf::st_transform("epsg:4326")
grid
#> Simple feature collection with 6 features and 20 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -95.99963 ymin: 22.99903 xmax: -95.99962 ymax: 22.99903
#> Geodetic CRS:  WGS 84
#>                Event_ID irwinID Incid_Name Incid_Type Map_ID Map_Prog
#> 0 CO3922010528720020608    <NA>     HAYMAN   Wildfire  12585     MTBS
#> 1 CO3922010528720020608    <NA>     HAYMAN   Wildfire  12585     MTBS
#> 2 CO3922010528720020608    <NA>     HAYMAN   Wildfire  12585     MTBS
#> 3 CO3922010528720020608    <NA>     HAYMAN   Wildfire  12585     MTBS
#> 4 CO3922010528720020608    <NA>     HAYMAN   Wildfire  12585     MTBS
#> 5 CO3922010528720020608    <NA>     HAYMAN   Wildfire  12585     MTBS
#>   Asmnt_Type BurnBndAc    Ig_Date          Pre_ID         Post_ID Perim_ID
#> 0   Extended    128726 2002/06/08 503303320010824 503303320030814     <NA>
#> 1   Extended       265 2002/06/08 503303320010824 503303320030814     <NA>
#> 2   Extended       322 2002/06/08 503303320010824 503303320030814     <NA>
#> 3   Extended        29 2002/06/08 503303320010824 503303320030814     <NA>
#> 4   Extended        50 2002/06/08 503303320010824 503303320030814     <NA>
#> 5   Extended        25 2002/06/08 503303320010824 503303320030814     <NA>
#>   dNBR_offst dNBR_stdDv NoData_T IncGreen_T Low_T Mod_T High_T Comment
#> 0         29      -9999     -970       -150   140   211    350    <NA>
#> 1         29      -9999     -970       -150   140   211    350    <NA>
#> 2         29      -9999     -970       -150   140   211    350    <NA>
#> 3         29      -9999     -970       -150   140   211    350    <NA>
#> 4         29      -9999     -970       -150   140   211    350    <NA>
#> 5         29      -9999     -970       -150   140   211    350    <NA>
#>                     geometry
#> 0 POINT (-95.99963 22.99903)
#> 1 POINT (-95.99962 22.99903)
#> 2 POINT (-95.99962 22.99903)
#> 3 POINT (-95.99962 22.99903)
#> 4 POINT (-95.99963 22.99903)
#> 5 POINT (-95.99962 22.99903)

You could visually overlay your data with e.g. plot(p); plot(grid, add = TRUE) but this won't really work with the subset provided since the extents do not overlap. grid (at least the first objects you provided) are located somehwere in the Gulf of Mexico, but this maybe is not an issue with your original dataset if you're covering a large area.

Upvotes: 0

Related Questions