search

ropensci / FedData

Functions to Automate Downloading Geospatial Data Available from Several Federated Data Sources
https://docs.ropensci.org/FedData
Other
99 stars 22 forks source link

Help Requested: Cropping NLCD_Data RasterLayer to Shapefiles #57

Closed fishcakebaker closed 4 years ago

fishcakebaker commented 4 years ago

Hi,

First, thank you for this package! It's a really great collection. Also thank you for the updated solution for the issue listed here, I had exactly the same problem and your solutions worked.

I'm just having a bit (massive) mental block trying to understand how the get_nlcd function is interpreting my shape file.

I'm using US Census 2017 data, specifically the LA County shape file as extracted in my code snippet below.

install.packages("devtools") devtools::install_github("ropensci/FedData") devtools::install_github("bocinsky/paleocar") library(FedData) library(sf) shp <- read_sf("tl_2017_us_county/tl_2017_us_county.shp") sub_shp <- subset(shp, GEOID == "06037") nlcd_data <- get_nlcd(sub_shp[18], year = 2001, label = "06037", dataset = "Land_Cover") plot(nlcd_data)

Here the GEOID is the 6 digit FIPS number (state 2 + county 4). This code outputs fine, with the image below. I pull sub_shp[18] as this is the geometry Column, though it is not actually required in get_nlcd().

NLCD_LACounty

My issue is that I want the resultant RasterLayer to be cropped to the profile/shape/outline of the Shapefiles, so I can then just carry out basic pixel counting from the legend entry.

I'm pretty new to R, as in only one or two weeks in. So I could just be missing something obvious.

I suspected that it would just need me to use the raster::crop function, e.g.

rtest1 <- raster::crop(nlcd_data,sub_shp[18])

However this outputs the error: Error in .local(x, y, ...) : extents do not overlap

From my Google, Crop seemed like the way to go. But am I missing something? Perhaps there's a transformation to map the NLCD layer across the shape file to crop it.

Sorry if I've not been clear enough! It's been a long week.

bocinsky commented 4 years ago

Hey there. I’m not really able to help too much today… constant zoom meetings… but I think the problem is likely that the coordinate reference systems aren’t the same between the NLCD output and your study area. You need to transform your shape file to the right CRS. Try:

rtest1 <- raster::crop(nlcd_data, sf::st_transform( sub_shp[18], raster::projection(nlcd_data)))

Let me know if that works for you.

On May 15, 2020, at 12:26 PM, Reece notifications@github.com wrote:

Hi,

First, thank you for this package! It's a really great collection. Also thank you for the updated solution for the issue listed here https://github.com/ropensci/FedData/issues/46#issue-479779287, I had exactly the same problem and your solutions worked.

I'm just having a bit (massive) mental block trying to understand how the get_nlcd function is interpreting my shape file.

I'm using US Census 2017 data, https://www.census.gov/geographies/mapping-files/time-series/geo/tiger-line-file.2017.html specifically the LA County shape file as extracted in my code snippet below.

install.packages("devtools") devtools::install_github("ropensci/FedData") devtools::install_github("bocinsky/paleocar") library(FedData) library(sf) shp <- read_sf("tl_2017_us_county/tl_2017_us_county.shp") sub_shp <- subset(shp, GEOID == "06037") nlcd_data <- get_nlcd(sub_shp[18], year = 2001, label = "06037", dataset = "Land_Cover") plot(nlcd_data)

Here the GEOID is the 6 digit FIPS number (state 2 + county 4). This code outputs fine, with the image below. I pull sub_shp[18] as this is the geometry Column, though it is not actually required in get_nlcd().

https://user-images.githubusercontent.com/47860255/82082843-a008e680-96e0-11ea-93d8-943702021c0c.png My issue is that I want the resultant RasterLayer to be cropped to the profile/shape/outline of the Shapefiles, so I can then just carry out basic pixel counting from the legend entry.

I'm pretty new to R, as in only one or two weeks in. So I could just be missing something obvious.

I suspected that it would just need me to use the raster::crop function, e.g.

rtest1 <- raster::crop(nlcd_data,sub_shp[18])

However this outputs the error: Error in .local(x, y, ...) : extents do not overlap

From my Google, Crop seemed like the way to go. But am I missing something? Perhaps there's a transformation to map the NLCD layer across the shape file to crop it.

Sorry if I've not been clear enough! It's been a long week.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/ropensci/FedData/issues/57, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB7SSM7TOMOIXLA2X4FDAT3RRWCORANCNFSM4NBZIQRQ.

fishcakebaker commented 4 years ago

Thanks for the reply! I wasn't expecting one so quick to be honest, I know the feeling with Zoom meetings.

I think it is the CRS issue, but the code snippet above doesn't change the image. I've posted my total code below, with outputs. 

devtools::install_github("ropensci/FedData")
devtools::install_github("bocinsky/paleocar")

library(FedData)
library(sf)

shp <- read_sf("tl_2017_us_county/tl_2017_us_county.shp")

sub_shp <- subset(shp, GEOID == "06037")

nlcd_data <- get_nlcd(sub_shp, 
                      year = 2016,
                      label = "06037",
                      dataset = "Land_Cover")

plot(nlcd_data)

rtest1 <- raster::crop(nlcd_data, sf::st_transform(sub_shp, raster::projection(nlcd_data)))
plot(rtest1)

Original Output image

New Output image

I was curious to the flow of transformation. Which just proved that the transformation is working, but I've pasted it below for good measure.

st_crs(sub_shp)

  User input: 4269 
  wkt:
GEOGCS["NAD83",
    DATUM["North_American_Datum_1983",
        SPHEROID["GRS 1980",6378137,298.257222101,
            AUTHORITY["EPSG","7019"]],
        TOWGS84[0,0,0,0,0,0,0],
        AUTHORITY["EPSG","6269"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4269"]]```

`st_crs(st_transform(sub_shp,st_crs(nlcd_data))`
```Coordinate Reference System:
  User input: +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs 
  wkt:
PROJCS["unnamed",
    GEOGCS["unnamed ellipse",
        DATUM["unknown",
            SPHEROID["unnamed",6378137,0],
            EXTENSION["PROJ4_GRIDS","@null"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Mercator_2SP"],
    PARAMETER["standard_parallel_1",0],
    PARAMETER["central_meridian",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]

Which is now the same CRS as the nlcd_data

st_crs(nlcd_data)

  User input: +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs 
  wkt:
PROJCS["unnamed",
    GEOGCS["unnamed ellipse",
        DATUM["unknown",
            SPHEROID["unnamed",6378137,0],
            EXTENSION["PROJ4_GRIDS","@null"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Mercator_2SP"],
    PARAMETER["standard_parallel_1",0],
    PARAMETER["central_meridian",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]

So the transformation works, but it doesn't affect the cropping. When I look at rtest1 and nlcd_data they also look identical.

bocinsky commented 4 years ago

Ah, cropping a raster just crops to a bounding box. You are probably looking for the raster::mask function, which will set pixels outside of the polygon to NA (making them transparent when plotting).

Cheers! KB

On May 15, 2020, at 4:52 PM, Reece notifications@github.com wrote:

Thanks for the reply! I wasn't expecting one so quick to be honest, I know the feeling with Zoom meetings.

I think it is the CRS issue, but the code snippet above doesn't change the image. I've posted my total code below, with outputs.

devtools::install_github("ropensci/FedData") devtools::install_github("bocinsky/paleocar")

library(FedData) library(sf)

shp <- read_sf("tl_2017_us_county/tl_2017_us_county.shp")

sub_shp <- subset(shp, GEOID == "06037")

nlcd_data <- get_nlcd(sub_shp, year = 2016, label = "06037", dataset = "Land_Cover")

plot(nlcd_data)

rtest1 <- raster::crop(nlcd_data, sf::st_transform(sub_shp, raster::projection(nlcd_data))) plot(rtest1) Original Output https://user-images.githubusercontent.com/47860255/82101684-16b7db00-9705-11ea-920c-83cbc949f6fd.png New Output https://user-images.githubusercontent.com/47860255/82101692-1c152580-9705-11ea-92ff-59e3a35a9d5e.png I was curious to the flow of transformation. Which just proved that the transformation is working, but I've pasted it below for good measure.

st_crs(sub_shp)

User input: 4269 wkt: GEOGCS["NAD83", DATUM["North_American_Datum_1983", SPHEROID["GRS 1980",6378137,298.257222101, AUTHORITY["EPSG","7019"]], TOWGS84[0,0,0,0,0,0,0], AUTHORITY["EPSG","6269"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4269"]]```

st_crs(st_transform(sub_shp,st_crs(nlcd_data))


  User input: +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs 
  wkt:
PROJCS["unnamed",
    GEOGCS["unnamed ellipse",
        DATUM["unknown",
            SPHEROID["unnamed",6378137,0],
            EXTENSION["PROJ4_GRIDS","@null"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Mercator_2SP"],
    PARAMETER["standard_parallel_1",0],
    PARAMETER["central_meridian",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
Which is now the same CRS as the nlcd_data

st_crs(nlcd_data)

  User input: +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs 
  wkt:
PROJCS["unnamed",
    GEOGCS["unnamed ellipse",
        DATUM["unknown",
            SPHEROID["unnamed",6378137,0],
            EXTENSION["PROJ4_GRIDS","@null"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Mercator_2SP"],
    PARAMETER["standard_parallel_1",0],
    PARAMETER["central_meridian",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
So the transformation works, but it doesn't affect the cropping. When I look at rtest1 and nlcd_data they also look identical.

—
You are receiving this because you commented.
Reply to this email directly, view it on GitHub <https://github.com/ropensci/FedData/issues/57#issuecomment-629538889>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AB7SSM75KSHOL3I5W46NM4DRRXBTTANCNFSM4NBZIQRQ>.
fishcakebaker commented 4 years ago

You're right, and it works! Just incase anyone else is looking at this later, below is the updated line and output.

Thank you for the virtual hand-holding through the problem. I'll try to pass it on one day. Hope you have a good weekend!

rtest1 <- raster::mask(nlcd_data, sf::st_transform(sub_shp, raster::projection(nlcd_data)))

image

bocinsky commented 4 years ago

Super! Happy mapping!