How To Change The Projection Of A Raster In R
Lesson half-dozen. GIS in R: How to Reproject Vector Information in Different Coordinate Reference Systems (crs) in R
Learning Objectives
Afterward completing this tutorial, you will be able to:
- Describe atleast 2 reasons that a information provider may chose to store a dataset in a item
CRS. - Reproject a vector dataset to another
CRSinR. - Identify the
CRSof a spatial dataset inR.
What You lot Demand
You will need a figurer with internet access to complete this lesson and the data for calendar week 4 of the grade.
Download Week 4 Information (~500 MB)
Working with Spatial Data from Unlike Sources
Yous often demand to gather spatial datasets for from dissimilar sources and/or data that embrace unlike spatial extents. Spatial data from different sources and that embrace different extents are frequently in different Coordinate Reference Systems (CRS).
Some reasons for data beingness in dissimilar CRSsouthward include:
- The data are stored in a detail
CRSconvention used by the data provider which might be a federal agency, or a land planning office. - The data are stored in a particular
CRSthat is customized to a region. For instance, many states prefer to use a Land Plane projection customized for that country.
In this tutorial you will larn how to place and manage spatial data in different projections. You will learn how to reproject the data then that they are in the aforementioned projection to support plotting / mapping. Note that these skills are besides required for any geoprocessing / spatial analysis. Data need to be in the same CRS to ensure accurate results.
You volition use the rgdal and raster libraries in this tutorial.
# load spatial data packages library ( rgdal ) library ( raster ) library ( rgeos ) options ( stringsAsFactors = FALSE ) # fix working directory to information binder # setwd("pathToDirHere") Import United states of america Boundaries - Census Information
There are many good sources of boundary base of operations layers that you lot tin use to create a basemap. Some R packages even have these base of operations layers congenital in to support quick and efficient mapping. In this tutorial, you volition utilize boundary layers for the Usa, provided by the The states Demography Agency.
It is useful to accept shapefiles to piece of work with considering yous can add additional attributes to them if need exist - for projection specific mapping.
Read US Purlieus File
You will use the readOGR() function to import the /united states-purlieus-layers/United states of america-Land-Boundaries-Demography-2014 layer into R. This layer contains the boundaries of all continental states in the U.S. Delight note that these data have been modified and reprojected from the original data downloaded from the Census website to back up the learning goals of this tutorial.
# Import the shapefile information into R state_boundary_us <- readOGR ( "information/week-04/the states-boundary-layers" , "The states-State-Boundaries-Census-2014" ) ## OGR data source with driver: ESRI Shapefile ## Source: "/root/earth-analytics/data/week-04/u.s.a.-purlieus-layers", layer: "Us-Country-Boundaries-Demography-2014" ## with 58 features ## It has x fields ## Integer64 fields read as strings: ALAND AWATER # view data construction class ( state_boundary_us ) ## [1] "SpatialPolygonsDataFrame" ## attr(,"package") ## [ane] "sp" Note: The Z-dimension alert is normal. The readOGR() function doesn't import z (vertical dimension or height) data by default. This is considering not all shapefiles contain z dimension data. More on readOGR
Next, let'due south plot the U.Southward. states data.
# view cavalcade names plot ( state_boundary_us , main = "Map of Continental U.s. State Boundaries\n Usa Census Bureau Information" ) 
U.S. Boundary Layer
You can add together a boundary layer of the United States to your map - to make information technology look nicer. Yous volition import data/week-04/usa-purlieus-layers/US-Boundary-Dissolved-States. If you specify a thicker line width using lwd = iv for the border layer, it will make your map popular!
# Read the .csv file country_boundary_us <- readOGR ( "data/week-04/usa-boundary-layers" , "US-Boundary-Dissolved-States" ) ## OGR data source with commuter: ESRI Shapefile ## Source: "/root/earth-analytics/data/week-04/u.s.-boundary-layers", layer: "U.s.-Boundary-Dissolved-States" ## with 1 features ## It has 9 fields ## Integer64 fields read as strings: ALAND AWATER # look at the data structure grade ( country_boundary_us ) ## [1] "SpatialPolygonsDataFrame" ## attr(,"bundle") ## [1] "sp" # view cavalcade names plot ( state_boundary_us , main = "Map of Continental US State Boundaries\n US Census Bureau Data" , border = "gray40" ) # view column names plot ( country_boundary_us , lwd = 4 , edge = "gray18" , add = TRUE ) 
Side by side, let's add the location of your study area sites. As you are adding these layers, take annotation of the class of each object. You will use AOI to correspond "Area of Interest" in your data.
# Import a polygon shapefile sjer_aoi <- readOGR ( "data/week-04/california/SJER/vector_data" , "SJER_crop" ) ## OGR data source with driver: ESRI Shapefile ## Source: "/root/earth-analytics/data/calendar week-04/california/SJER/vector_data", layer: "SJER_crop" ## with i features ## It has one fields class ( sjer_aoi ) ## [1] "SpatialPolygonsDataFrame" ## attr(,"package") ## [1] "sp" # plot indicate - looks ok? plot ( sjer_aoi , pch = 19 , col = "purple" , chief = "San Joachin Experimental Range AOI" ) 
Your SJER AOI layer plots nicely. Let'due south next add it equally a layer on top of the U.Southward. states and boundary layers in your basemap plot.
# plot state boundaries plot ( state_boundary_us , main = "Map of Continental US Country Boundaries \north with SJER AOI" , border = "gray40" ) # add United states border outline plot ( country_boundary_us , lwd = 4 , border = "gray18" , add = TRUE ) # add AOI boundary plot ( sjer_aoi , pch = 19 , col = "majestic" , add together = TRUE ) 
What do you detect about the resultant plot? Do y'all see the AOI boundary in the California area? Is something wrong?
Let's bank check out the CRS (crs()) of both datasets to see if you can identify any problems that might cause the point location to non plot properly on top of your U.S. boundary layers.
# view CRS of your site data crs ( sjer_aoi ) ## CRS arguments: ## +proj=utm +zone=eleven +datum=WGS84 +units=one thousand +no_defs +ellps=WGS84 ## +towgs84=0,0,0 # view crs of census information crs ( state_boundary_us ) ## CRS arguments: ## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 crs ( country_boundary_us ) ## CRS arguments: ## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 It looks like your data are in different CRS. You tin can tell this by looking at the CRS strings in proj4 format.
Understanding CRS in proj4 Format
The CRS for your information are given to the states past R in proj4 format. Let'due south break down the pieces of proj4 string. The cord contains all of the individual CRS elements that R or some other GIS might need. Each element is specified with a + sign, similar to how a .csv file is delimited or broken upwards by a ,. After each + you see the CRS chemical element being defined. For instance projection (proj=) and datum (datum=).
UTM proj4 String
Your project string for sjer_aoi specifies the UTM projection as follows:
+proj=utm +zone=18 +datum=WGS84 +units=thou +no_defs +ellps=WGS84 +towgs84=0,0,0
- proj=utm: the projection is UTM, UTM has several zones
- zone=18: the zone is 18
- datum=WGS84: the datum WGS84 (the datum refers to the 0,0 reference for the coordinate arrangement used in the projection)
- units=g: the units for the coordinates are in METERS
- ellps=WGS84: the ellipsoid (how the Earth'southward roundness is calculated) for the data is WGS84
Note that the zone is unique to the UTM project. Not all CRS will accept a zone.
Geographic (lat / long) proj4 String
Your project string for state_boundary_us and country_boundary_us specifies the lat/long projection as follows:
+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
- proj=longlat: the data are in a geographic (latitude and longitude) coordinate system
- datum=WGS84: the datum WGS84 (the datum refers to the 0,0 reference for the coordinate organization used in the project)
- ellps=WGS84: the ellipsoid (how the earth's roundness is calculated) is WGS84
Note that at that place are no specified units above. This is considering this geographic coordinate reference system is in latitude and longitude which is most often recorded in Decimal Degrees.
Information tip: the last portion of each proj4 string is +towgs84=0,0,0 . This is a conversion factor that is used if a datum conversion is required. You will non deal with datums in this tutorial serial.
CRS units - View Object Extent
Adjacent, allow'due south view the extent or spatial coverage for the sjer_aoi spatial object compared to the state_boundary_us object.
# extent & crs for AOI extent ( sjer_aoi ) ## form : Extent ## xmin : 254570.6 ## xmax : 258867.4 ## ymin : 4107303 ## ymax : 4112362 crs ( sjer_aoi ) ## CRS arguments: ## +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 ## +towgs84=0,0,0 # extent & crs for object in geographic extent ( state_boundary_us ) ## class : Extent ## xmin : -124.7258 ## xmax : -66.94989 ## ymin : 24.49813 ## ymax : 49.38436 crs ( state_boundary_us ) ## CRS arguments: ## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 Note the divergence in the units for each object. The extent for state_boundary_us is in latitude and longitude which yields smaller numbers representing decimal degree units. Your AOI boundary indicate is in UTM, represented in meters.
Proj4 & CRS Resources
- More data on the proj4 format
- A fairly comprehensive list of CRS past format
- To view a list of datum conversion factors blazon:
projInfo(type = "datum")into theRpanel.
Reproject Vector Information
Now yous know your data are in dissimilar CRS. To address this, yous have to modify or reproject the data so they are all in the aforementioned CRS. You can apply spTransform() function to reproject your data. When yous reproject the data, you specify the CRS that you wish to transform your data to. This CRS contains the datum, units and other information that R needs to reproject your information.
The spTransform() function requires two inputs:
- the proper noun of the object that you wish to transform
- the
CRSthat you wish to transform that object also. In this example you can employ thecrs()of thestate_boundary_usobject equally follows:crs(state_boundary_us)
Data Tip: spTransform() will simply piece of work if your original spatial object has a CRS assigned to it AND if that CRS is the correct CRS!
Adjacent, let's reproject your betoken layer into the geographic - breadth and longitude WGS84 coordinate reference organisation (CRS).
# reproject data sjer_aoi_WGS84 <- spTransform ( sjer_aoi , crs ( state_boundary_us )) # what is the CRS of the new object crs ( sjer_aoi_WGS84 ) ## CRS arguments: ## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 # does the extent look like decimal degrees? extent ( sjer_aoi_WGS84 ) ## class : Extent ## xmin : -119.7626 ## xmax : -119.7127 ## ymin : 37.0799 ## ymax : 37.12657 One time your data are reprojected, you can endeavour to plot again.
# plot state boundaries plot ( state_boundary_us , main = "Map of Continental The states State Boundaries\due north With SJER AOI" , edge = "gray40" ) # add together US border outline plot ( country_boundary_us , lwd = 4 , border = "gray18" , add = TRUE ) # add AOI plot ( sjer_aoi_WGS84 , pch = 19 , col = "purple" , add together = Truthful ) 
But at present, the AOI is a polygon and it'southward too small to encounter on the map. Let'southward convert the polygon to a polygon CENTROID and plot yet over again.
# get coordinate center of the polygon aoi_centroid <- coordinates ( sjer_aoi_WGS84 ) # plot country boundaries plot ( state_boundary_us , main = "Map of Continental US Land Boundaries\n With SJER AOI" , edge = "gray40" ) # add US border outline plot ( country_boundary_us , lwd = 4 , edge = "gray18" , add = Truthful ) # add betoken location of the centroid to the plot points ( aoi_centroid , pch = 8 , col = "magenta" , cex = 1.5 ) 
Reprojecting your information ensured that things lined up on your map! Information technology will besides allow us to perform any required geoprocessing (spatial calculations / transformations) on your data.
Exam Your Skills: Crop, Peproject, Plot Information
Create a map of your SJER study area as follows:
- Import the
madera-county-roads/tl_2013_06039_roads.shplayer located in your week4 data download. - Create a map that shows the roads layer, study site locations and the
sjer_aoiboundary. - Add a title to your plot.
- Add a fable to your plot that shows both the roads and the plot locations.
- Plot the roads by road blazon and add together each type to the fable. HINT: use the metadata included in your information download to figure out what each type of road represents ("C", "S", etc.). Utilise the homework lesson on custom legends to assistance build the legend.
- BONUS: Plot the plots by type - arrange the symbology of the plot locations (cull a symbol using pch for each type and adjust the color of the points).
- Practice your best to make the map await nice!
IMPORTANT: exist sure that all of the data are inside the same EXTENT and crs of the sjer_aoi layer. This means that you may take to CROP and reproject your data prior to plotting it!
Your map should look something like the map beneath. You should of grade use the actual roads types that y'all discover in the metadata rather than "Route type ane, etc"
NOTE: this is also a plot you will submit every bit a part of your homework this week!
Source: https://www.earthdatascience.org/courses/earth-analytics/spatial-data-r/reproject-vector-data/
Posted by: robbinssencong.blogspot.com
0 Response to "How To Change The Projection Of A Raster In R"
Post a Comment