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

Question

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

Welcome! This sounds like a pretty basic task and should be solved within minutes. Unfortunately, right now, your code is not really reproducible because we do not know what `Robindf` looks like. Would you please run `head(Robindf, 10) |> dput()` and edit your question with the object returned? Have a look at [how to make a great R reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) as reference. — dimfalk, Oct 28 '22 at 06:39
Is your wildlife grid a raster or a vector dataset consisting of polygons, e.g. a shapefile?`dput()` on raster objects is rather difficult - here it would be best to upload the file externally (dropbox/drive) and provide a link. — dimfalk, Oct 28 '22 at 18:21
My wildlife grid is a vector dataset (that I created myself), and I'm trying to display it on a shapefile grid that I downloaded from a wildfire database. I am having trouble getting the crs of both to align, if that helps clarify things at all. Thanks again for reading this. — BiggestNoob, Oct 28 '22 at 18:26
Then probably it would be best if you also provided these datasets to be able to troubleshoot properly. — dimfalk, Oct 28 '22 at 18:28
Thank you, I've gone in and added some information about that shapefile. — BiggestNoob, Oct 28 '22 at 18:37

score 0 · Answer 1 · answered Oct 28 '22 at 19:15

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.

Thank you!! I am indeed working with points near Denver. I was incorrectly inputting the "crs = "epsg:4326" portion each time. Thank you for helping me sort that. Confusingly enough, the hayman file should cover an extent in the same area. I'm not sure why when reprojected it shows up with such odd coordinates, but oh well. This helps a lot. — BiggestNoob, Oct 28 '22 at 19:40
You're welcome! If you created `hayman` yourself in e.g. ArcMap, maybe there was some confusion in crs definitions from the beginning when you started digitizing? Make sure the projection of your layer equals the one of your project, otherwise coordinates would not align, afaik. — dimfalk, Oct 28 '22 at 20:20
Thank you again! I actually did not create it, I downloaded it from mtbs.gov and there are many different layers in their dataset, so I will investigate whether there is a mismatch between those projections. — BiggestNoob, Oct 29 '22 at 01:17

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

1 Answers1