L15: Rasters and Raster Operations 1

Reading Rasters, Attributes, Raster Subsetting, Vectorizing, Cropping

Bogdan G. Popescu
bogdan.popescu@johncabot.edu

John Cabot University

Today

Today’s agenda includes:

  • becoming familiar with others types of data in R:
  • raster layers
  • examining spatial and non-spatial properties of vectors

We will use sf, stars, and units packages.

Reminder: Spatial Data Classes in R

  • Spatial entries in GIS can be represented as a vector or as a raster

  • In the previous lecture, we examine vectors and geometric operations with vectors.

Raster

Rasters are matrices or an array representing a rectangular area on the surface of the earth.

They contain elements such as :

  • Origin (\(x_{max}\), \(y_{min}\)) and resolution (\(\delta_{x}\), \(\delta_{y}\))
  • Coordinate Reference System
  • Values
  • Dimensions (rows, column, layers)
  • Extent: \(x_{min}\), \(y_{min}\), \(x_{max}\), \(y_{max}\)

Raster

Raster with the Same Extent but Different Resolutions

Raster file formats

  • “Simple” rasters
    • GeoTIFF (.tif)
    • Erdas Imagine Image (.img)
  • “Complex” rasters (>3D and / or metadata)
    • HDF (.hdf, .he5)
    • NetCDF (.nc)

Using R packages for rasters

There is one important package that we will work with: stars

terra is also another important package developed by the same authors

Both contain classes, and extensive functions for working with rasters in R

We will mostly work with stars

The stars package

Unlike terra and raster, stars is very well integrated with sf

sf is the package that we have used for vector analysis.

Example of Raster: Luminosity

A typical example of raster is satellite luminosity

Within economics, researchers use satellite luminosity as a proxy for economic activity

Luminosity in 2013

Working with Luminosity

You can download satellite luminosity by going to https://www.ncei.noaa.gov/products/dmsp-operational-linescan-system


Working with Luminosity

The next step is to scroll down to products


Working with Luminosity

The next step is to scroll down to products


Working with Luminosity

You can download these files by clicking on the links


Working with Luminosity

Note that these are large files: Every zip file will range from 355 MB to 411.1 MB

The files are cloud-free composites

In cases where two satellites were collecting data - two composites were produced.

Working with Luminosity


Working with Luminosity

Note that these are large files: Every zip file will range from 355 MB to 411.1 MB

The files are cloud-free composites

In cases where two satellites were collecting data - two composites were produced.

The products are 30 arc second grids, spanning -180 to 180 degrees longitude and -65 to 75 degrees latitude.

Each composite set is named with the satellite and the year (F121995 is from DMSP satellite number F12 for the year 1995)

Working with Luminosity

Three types of images are available:

  • F1?YYYY_v4b_cf_cvg.tif: Cloud-free coverages tally the total number of observations that went into each 30 arc second grid cell
  • F1?YYYY_v4b_avg_vis.tif Raw avg_vis contains the average of the visible band digital number values with no further filtering.
  • F1?YYYY_v4b_stable_lights.avg_vis.tif The cleaned up avg_vis contains the lights from cities, towns, and other sites with persistent lighting, including gas flares.

Working with Luminosity

The DMSP/OLS resolution is approximately 1000m

It comprises six different DMSP satellites:

  • F10 (1992–1994)
  • F12 (1994–1999)
  • F14 (1997–2003)
  • F15 (2000–2007)
  • F16 (2004–2009)
  • F18 (2010–2013)

Working with Luminosity

After downloading, your folder should look like below:


Working with Luminosity

You download all the relevant files for this exercise below:

Working with Luminosity

These are the steps to have the data ready for plotting in ggplot

#Step1: Read country borders
sf_use_s2(FALSE)
borders <- st_read(dsn="./data/boundaries/all-country-boundaries.shp", quiet = TRUE)
borders<-st_simplify(borders, dTolerance=0.1)

This what the shape file looks like:

ggplot()+
  geom_sf(data=borders)

Working with Luminosity

These are the steps to have the data ready for plotting in ggplot

#Step2: Creating a folder called "F101992"
target_dir <- "./data/luminosity/F101992"
#Step3: Create the target directory if it doesn't exist
if (!dir.exists(target_dir)) {
  dir.create(target_dir)
}

Working with Luminosity

This is what our folder looked like before running the code.


Working with Luminosity

These are the steps to have the data ready for plotting in ggplot

#Step2: Creating a folder called "F101992"
target_dir <- "./data/luminosity/F101992"
#Step3: Create the target directory if it doesn't exist
if (!dir.exists(target_dir)) {
  dir.create(target_dir)
}

Working with Luminosity

This is what our folder looks like after running the code.


Working with Luminosity

#Step4: Unzip the file
untar("./data/luminosity/F101992.v4.tar", exdir = target_dir)

Working with Luminosity

This is what our folder looks like before running the code.


Working with Luminosity

This is what our folder looks like after running the code.


Working with Luminosity

#Step5: Unzip the file further
R.utils::gunzip("./data/luminosity/F101992/F101992.v4b_web.avg_vis.tif.gz", remove = FALSE, overwrite=TRUE)

Working with Luminosity

This is what our folder looks like before running the code.


Working with Luminosity

This is what our folder looks like after running the code.


Working with Luminosity

#Step6: Read the raster
f1992<-read_stars("./data/luminosity/F101992/F101992.v4b_web.avg_vis.tif")
#Step7: Make sure the two files have the same CRS
st_crs(f1992)<-st_crs(borders)
plot(f1992)

Working with Luminosity

Note: do not try to plot it yet with ggplot.

The raster is too detailed.

It will freeze your computer.

Working with Luminosity

This is how we can reduce the resolution of a raster

#Step8: Reduce the resolution of the raster: Higher dx/day values mean lower resolution
f1992b<-st_warp(f1992, st_as_stars(st_bbox(f1992), dx = 0.7, dy = 0.7))
#Step8: Reduce the resolution of the raster: Higher dx/day values mean lower resolution
f1992b_reduced<-st_warp(f1992, st_as_stars(st_bbox(f1992), dx = 2.5, dy = 2.5))

Working with Luminosity

Let us now compare the two

Show the code 
ggplot()+
  geom_stars(data=f1992b)+
  scale_fill_gradientn(colours=c("black","white"))

Show the code 
ggplot()+
  geom_stars(data=f1992b_reduced)+
  scale_fill_gradientn(colours=c("black","white"))

Working with Luminosity

We can finally add a more friendly name to the raster

#Step9: Rename the raster
names(f1992b)<-"lights_1992"

Working with Luminosity

Just reviewing the steps:

#Step1: Read country borders
borders <- st_read(dsn="./data/boundaries/all-country-boundaries.shp", quiet = TRUE)
#Step2: Creating a folder called "F101992"
target_dir <- "./data/luminosity/F101992"
#Step3: Create the target directory if it doesn't exist
if (!dir.exists(target_dir)) {
  dir.create(target_dir)
}
#Step4: Unzip the file
untar("./data/luminosity/F101992.v4.tar", exdir = target_dir)
#Step5: Unzip the file further
R.utils::gunzip("./data/luminosity/F101992/F101992.v4b_web.avg_vis.tif.gz", remove = FALSE, overwrite=TRUE)
#Step6: Read the raster
f1992<-read_stars("./data/luminosity/F101992/F101992.v4b_web.avg_vis.tif")
#Step7: Make sure the two files have the same CRS
st_crs(f1992)<-st_crs(borders)
#Step8: Reduce the resolution of the raster: Higher dx/day values mean lower resolution
f1992b<-st_warp(f1992, st_as_stars(st_bbox(f1992), dx = 0.7, dy = 0.7))
#Step9: Rename the raster
names(f1992b)<-"lights_1992"

fig<-ggplot()+
  geom_stars(data=f1992b)+
  scale_fill_gradientn(colours=c("black","white"))+
  geom_sf(data=borders, linewidth = 0.1, fill = NA, color = "white", alpha=0.5)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Luminosity in 1992

Luminosity in 1992 in Europe

Show the code 
norway<-subset(borders, SOVEREIGNT == "Norway")
max_lat_y_eu<-st_bbox(norway)["ymax"]

greece<-subset(borders, SOVEREIGNT == "Greece")
min_lat_y_eu<-st_bbox(greece)["ymin"]

ukraine<-subset(borders, SOVEREIGNT == "Ukraine")
max_lon_x_eu<-st_bbox(ukraine)["xmax"]

portugal<-subset(borders, SOVEREIGNT == "Portugal")
min_lon_x_eu<-st_bbox(portugal)["xmin"]

ggplot()+
  geom_stars(data=f1992b)+
  scale_fill_gradientn(colours=c("black","white"))+
  geom_sf(data=borders, linewidth = 0.1, fill = NA, color = "white", alpha=0.5)+
  coord_sf(xlim = c(min_lon_x_eu-3, max_lon_x_eu+3), ylim = c(min_lat_y_eu-1, max_lat_y_eu-6))+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Raster Values and Properties

The print method gives a summary of the raster’s properties.

print(f1992b)
stars object with 2 dimensions and 1 attribute
attribute(s):
             Min. 1st Qu. Median     Mean 3rd Qu. Max. NA's
lights_1992     0       2      3 3.384738       4   63  515
dimension(s):
  from  to offset delta refsys x/y
x    1 515   -180   0.7 WGS 84 [x]
y    1 201     75  -0.7 WGS 84 [y]

The class function allows us to see the type of object that we’re dealing with

class(f1992b)
[1] "stars"

Raster Values and Properties

Fundamentally, a stars object is a collection of matrix or array objects.

This is along with additional properties of the dimensions, such as dimension names, Coordinate Reference Systems (CRS).

We can display the structure of a stars object using str.

str(f1992b)
List of 1
 $ lights_1992: num [1:515, 1:201] 0 0 0 0 0 0 0 0 0 0 ...
 - attr(*, "dimensions")=List of 2
  ..$ x:List of 7
  .. ..$ from  : num 1
  .. ..$ to    : num 515
  .. ..$ offset: num -180
  .. ..$ delta : num 0.7
  .. ..$ refsys:List of 2
  .. .. ..$ input: chr "WGS 84"
  .. .. ..$ wkt  : chr "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223"| __truncated__
  .. .. ..- attr(*, "class")= chr "crs"
  .. ..$ point : logi NA
  .. ..$ values: NULL
  .. ..- attr(*, "class")= chr "dimension"
  ..$ y:List of 7
  .. ..$ from  : num 1
  .. ..$ to    : num 201
  .. ..$ offset: num 75
  .. ..$ delta : num -0.7
  .. ..$ refsys:List of 2
  .. .. ..$ input: chr "WGS 84"
  .. .. ..$ wkt  : chr "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223"| __truncated__
  .. .. ..- attr(*, "class")= chr "crs"
  .. ..$ point : logi NA
  .. ..$ values: NULL
  .. ..- attr(*, "class")= chr "dimension"
  ..- attr(*, "raster")=List of 4
  .. ..$ affine     : num [1:2] 0 0
  .. ..$ dimensions : chr [1:2] "x" "y"
  .. ..$ curvilinear: logi FALSE
  .. ..$ blocksizes : NULL
  .. ..- attr(*, "class")= chr "stars_raster"
  ..- attr(*, "class")= chr "dimensions"
 - attr(*, "class")= chr "stars"

Raster Values and Properties

From this, we see that we have the f1992b object of length 1.

The 2nd line in the str printout specifies the name and contents of the first (and only) item in the list, namely its type, dimensions, and a sample of the first few values.

In our case, the first item is lights_1992, and it is a matrix with 515 rows and 201 columns

Raster Attributes and Values

To access some of the attributes of our objects, we can type:

names(f1992b)
[1] "lights_1992"

We can change names easily in the following wat:

names(f1992b)<-"lum_1992"
names(f1992b)
[1] "lum_1992"

Accessing Attributes by Name

We can access our raster’s attributes by using the list access methods

#f1992b$lights_1992
f1992b[["lum_1992"]][1:10, 1:10]
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    0    0    0    0    0    0    0   10    7     9
 [2,]    0    0    0    0    0    0    0    6    7     6
 [3,]    0    0    0    0    0    0    0    5    6     9
 [4,]    0    0    0    0    0    0    0    9    6     9
 [5,]    0    0    0    0    0    0    0    8    8     8
 [6,]    0    0    0    0    0    0    0    9    9     7
 [7,]    0    0    0    0    0    0    0    8    7     8
 [8,]    0    0    0    0    0    0    0    4   10     9
 [9,]    0    0    0    0    0    0    0    6    7     9
[10,]    0    0    0    0    0    0    0    0    8     9

This is how we demonstrate that our file is a matrix

#Using names
class(f1992b$lum_1992)
[1] "matrix" "array" 
#Using index
class(f1992b[[1]])
[1] "matrix" "array"

Accessing Attributes by Name

And that it has 512 rows and 201 columns.

#Using names
dim(f1992b$lum_1992)
  x   y 
515 201 
#Using index
dim(f1992b[[1]])
  x   y 
515 201

Subset operators in R

Syntax Objects Returns
x[i] vector, list Subset i
x[i, j] data.frame, matrix Subset i, j
x[i, j, k] array Subset i, j, k
x[[i]] vectors, lists Single element i
x$i data.frame, list Single column or element i
x@n S4 objects Slot n

Dimensions and spatial properties

We can use the following function to get at important raster properties

# Number of rows
nrow(f1992b)
  x 
515 
# Number of columns
ncol(f1992b)
  y 
201 
# All dimensions
dim(f1992b)
  x   y 
515 201 
# Extent of a raster
st_bbox(f1992b)
      xmin       ymin       xmax       ymax 
-180.00417  -65.69583  180.49583   75.00417

Dimensions and spatial properties

The other properties of rasters can be extracted in the following way:

st_dimensions(f1992b)$x$offset  ## x-axis origin
[1] -180.0042
st_dimensions(f1992b)$y$offset  ## x-axis origin
[1] 75.00417
st_dimensions(f1992b)$x$delta  ## x-axis resolution
[1] 0.7
st_dimensions(f1992b)$y$delta  ## y-axis resolution
[1] -0.7

Note that the resolution is separate for x and y.

If the (absolute) resolution is equal for both, raster pixels are square.

When the x and y resolutions are unequal, raster pixels are non-square rectangles.

Dimensions and spatial properties

We can find out the crs of our raster using the st_crs function

st_crs(f1992b)
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

Accessing Raster Values

As indicated, raster values can be extracted directly, as a matrix or as an array.

For example:

#Printing only the first 10 rows and the first 10 columns
f1992b$lum_1992[1:10, 1:10]
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    0    0    0    0    0    0    0   10    7     9
 [2,]    0    0    0    0    0    0    0    6    7     6
 [3,]    0    0    0    0    0    0    0    5    6     9
 [4,]    0    0    0    0    0    0    0    9    6     9
 [5,]    0    0    0    0    0    0    0    8    8     8
 [6,]    0    0    0    0    0    0    0    9    9     7
 [7,]    0    0    0    0    0    0    0    8    7     8
 [8,]    0    0    0    0    0    0    0    4   10     9
 [9,]    0    0    0    0    0    0    0    6    7     9
[10,]    0    0    0    0    0    0    0    0    8     9

Accessing Raster Values

A histogram will help us identify the raster values distribution

#Turning the object into a dataframe
df<-st_as_sf(f1992b)

Here is the difference between the original raster and the dataframe

#The raster
class(f1992b)
[1] "stars"
#The dataframe
class(df)
[1] "sf"         "data.frame"

Accessing Raster Values

A histogram will help us identify the raster values distribution

And here is how we can plot the raster vs. the dataframe

Show the code 
ggplot()+
  geom_stars(data=f1992b)+
  scale_fill_gradientn(colours=c("black","white"))

Show the code 
ggplot()+
  geom_sf(data=df, aes(fill=lum_1992),
          linewidth=0)+
  scale_fill_gradientn(colours=c("black","white"))+
  theme_bw()

Accessing Raster Values

And this is what the data looks like:

#Creating a histogram
ggplot(df, aes(x = lum_1992)) +
  geom_histogram(binwidth = 10, color = "white") +
  ggtitle("Histogram of lum_1992") +
  xlab("lum_1992") +
  ylab("Frequency")+
  theme_bw()

Accessing Raster Values

We thus see that most of values around the world are 0.

This is because 71 percent of the Earth’s surface is water-covered.

The important part is that by transforming the stars object into a dataframe (grid), we can perform the usual operations dedicated to dataframes.

For example: glimpse:

library(dplyr)
glimpse(df)
Rows: 103,000
Columns: 2
$ lum_1992 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ geometry <POLYGON [°]> POLYGON ((-180.0042 75.0041..., POLYGON ((-179.3042 7…

Accessing Raster Values

We can also calculate the mean, sd, min, and max

mean(df$lum_1992, na.rm=T)
[1] 3.384738
sd(df$lum_1992, na.rm=T)
[1] 3.161684
min(df$lum_1992, na.rm=T)
[1] 0
max(df$lum_1992, na.rm=T)
[1] 63

Mapping the World

This is what luminosity around the world looks like

Show the code 
p1<- ggplot()+
  geom_sf(data=df, aes(fill=lum_1992),
          linewidth=0)+
  scale_fill_gradientn(colours=c("black","white"))+
  geom_sf(data=borders, linewidth = 0.1, fill = NA, color = "white", alpha=0.5)+
  xlab("Longitude") + ylab("Latitude")+
  theme_bw()+
  coord_sf()
p1

The Histogram

This is what the data looks like:

Show the code 
p2<- ggplot(df, aes(lum_1992)) + 
  guides(color = "none") +
  geom_histogram(binwidth = 2, 
                 color = "white",
                 alpha = 0.5)+
  #Adding density and multiplying by a manually defined value
  #for visibility
  geom_density(aes(y = after_stat(density)* 80000), 
               fill=NA, color=alpha("blue", 0.5))+
  #Adding a secondary axis
  scale_y_continuous(name = "count",
                     sec.axis = 
                       sec_axis(~.x/80000, 
                                name = "density"))+
  theme_bw()
p2

Mapping Lights in Central Rome

Let’s focus on Italy.

You can download the Rome central area polygon

You can also download the main streets for the central area

Make sure you place them in the relevant folder.

Mapping Lights in Central Rome: Cropping Rasters

Let’s focus on Italy

#Step1: Reading the central polygon
relev_poly<-st_read(dsn="./data/rome_center/central_polygon.shp", quiet = T)
#Step2: Reading roads
gis_osm_roads <- st_read(dsn="./data/selection.gdb", layer="gis_osm_roads", quiet = T)
#Step3: Extracting the bix from the central polygon
ext2<-st_bbox(relev_poly)
#Step4: Ensuring the same crs for the polygon as for the raster
st_crs(relev_poly)<-st_crs(f1992)
#Step5: Cropping the raster
rc <- st_crop(f1992, relev_poly)

Mapping Lights in Central Rome: Cropping Rasters

This is how the original vs. the cropped raster compare

Show the code 
ggplot()+
  geom_stars(data=f1992b)+
  scale_fill_gradientn(colours=c("black","white"))

Show the code 
ggplot()+
  geom_stars(data=rc)+
  scale_fill_gradientn(colours=c("black","white"), na.value=NA)

Mapping Lights in Central Rome

We can finally add the roads in Rome to get more context.

#Step6: Mapping everything
ggplot()+
  geom_stars(data = rc)+
  geom_sf(data = relev_poly, color="white", fill=NA)+
  geom_sf(data = gis_osm_roads, color="grey50", linewidth=0.05)+
  theme_bw()+
  scale_fill_gradientn(colours=c("black","white"), na.value=NA)

Mapping Lights in Central Rome

Let’s focus on Italy

Mapping Lights in Central Rome

Every grid has a value:

#Turning the object into a dataframe
rc2<-st_as_sf(rc)
names(rc2)[1]<-"lum_1992"
#Step6: Mapping everything
ggplot()+
  geom_sf(data = rc2, aes(fill=lum_1992), linewidth=0)+
  geom_sf(data = relev_poly, color="white", fill=NA)+
  geom_sf(data = gis_osm_roads, color="grey50", linewidth=0.05)+
  theme_bw()+
  scale_fill_gradientn(colours=c("black","white"), na.value=NA)+
  geom_sf_text(data = rc2, 
  aes(label = lum_1992, geometry = geometry),
  color = ifelse(rc2$lum_1992 < 61, "white", "black"), size = 1.5)

Mapping Lights in Central Rome

Every grid has a value:

Mapping Lights in Central Rome

Let’s zoom in:

#Turning the object into a dataframe
rc2<-st_as_sf(rc)
names(rc2)[1]<-"lum_1992"
#Step6: Mapping everything
ggplot()+
  geom_sf(data = rc2, aes(fill=lum_1992), linewidth=0)+
  geom_sf(data = relev_poly, color="white", fill=NA)+
  geom_sf(data = gis_osm_roads, color="grey50", linewidth=0.05)+
  theme_bw()+
  scale_fill_gradientn(colours=c("black","white"), na.value=NA)+
  geom_sf_text(data = rc2, 
  aes(label = lum_1992, geometry = geometry),
  color = ifelse(rc2$lum_1992 < 61, "white", "black"),
  #aes(label = lum_1992, geometry = geometry, color = ifelse(lum_1992 < 61, "white", "black")),
  size = 4)+
  coord_sf(xlim = c(12.45, 12.60), 
           ylim = c(41.87, 41.79))

Mapping Lights in Central Rome

Let’s zoom in:

Mapping Lights in Central Rome

As you might imagine, we can use different colors for these values:

Mapping Lights in Central Rome

As you might imagine, we can use different colors for these values:

Interractive Maps with mapview

library(mapview)
mapview(rc2)

Replacing Values

We can easily replace the values inside the raster (grid) by using simple dataframe operations.

For example, let’s replace all the values that have 63 with NA:

rc3<-rc2
rc3$lum_1992[rc3$lum_1992==63]<-NA
pic<-ggplot()+
  geom_sf(data = rc3, aes(fill=lum_1992), linewidth=0)+
  geom_sf(data = relev_poly, color="white", fill=NA)+
  geom_sf(data = gis_osm_roads, color="grey50", linewidth=0.05)+
  theme_bw()+
  scale_fill_viridis_c(option = "magma",begin = 0.1)+
  #scale_fill_gradientn(colours=c("black","white"))+
  geom_sf_text(data = rc2, 
  aes(label = lum_1992, geometry = geometry),
  color = ifelse(rc3$lum_1992 < 61, "white", "black"),
  #aes(label = lum_1992, geometry = geometry, color = ifelse(lum_1992 < 61, "white", "black")),
  size = 4)+
  coord_sf(xlim = c(12.45, 12.60), 
           ylim = c(41.87, 41.79))
pic

Replacing Values

We can easily replace the values inside the raster (grid) by using simple dataframe operations.

For example, let’s replace all the values that have 63 with NA:

Replacing Values

We can easily replace the values inside the raster (grid) by using simple dataframe operations.

For example, let’s replace all the values that have 63 with NA:

Let us now make the NA grids transparent

rc3<-rc2
rc3$lum_1992[rc3$lum_1992==63]<-NA
pic<-ggplot()+
  geom_sf(data = rc3, aes(fill=lum_1992), linewidth=0)+
  geom_sf(data = relev_poly, color="white", fill=NA)+
  geom_sf(data = gis_osm_roads, color="grey50", linewidth=0.05)+
  theme_bw()+
  scale_fill_viridis_c(option = "magma",begin = 0.1,
                       na.value = "white")+
  #scale_fill_gradientn(colours=c("black","white"))+
  geom_sf_text(data = rc2, 
  aes(label = lum_1992, geometry = geometry),
  color = ifelse(rc3$lum_1992 < 61, "white", "black"),
  #aes(label = lum_1992, geometry = geometry, color = ifelse(lum_1992 < 61, "white", "black")),
  size = 4)+
  coord_sf(xlim = c(12.45, 12.60), 
           ylim = c(41.87, 41.79))
pic

Replacing Values

We can easily replace the values inside the raster (grid) by using simple dataframe operations.

For example, let’s replace all the values that have 63 with NA:

Let us now make the NA grids transparent

Writing raster to file

We can of course write the grid as a shapefile for we can turn it back to a raster.

rc3_raster = st_rasterize(rc3)

The raster

class(rc3_raster)
[1] "stars"

The grid

class(rc3)
[1] "sf"         "data.frame"

Writing raster to file

We can write the raster in the following way:

#Step1: Creating a folder called
target_dir <- "./data/output"
#Step2: Create the target directory if it doesn't exist
if (!dir.exists(target_dir)) {
  dir.create(target_dir)
}
#Step3: Writing the raster
write_stars(rc3_raster, './data/output/rc3_raster.tif')

Writing the grid to file

We can write the grid in the following way:

#Step1: Creating a folder called
target_dir <- "./data/output"
#Step2: Create the target directory if it doesn't exist
if (!dir.exists(target_dir)) {
  dir.create(target_dir)
}
#Step4: Writing the shape files
st_write(rc3, "./data/output/rc3_grid.shp", update = TRUE)
Updating layer `rc3_grid' to data source `./data/output/rc3_grid.shp' using driver `ESRI Shapefile'
Updating existing layer rc3_grid
Writing 535 features with 1 fields and geometry type Polygon.

Conclusion

We have covered a lot of operations related to rasters

We learned about:

  • different types of rasters
  • plotting rasters
  • dimensions of rasters
  • accessing raster values
  • working with luminosity
  • subsetting rasters
  • cropping rasters

Popescu (JCU): Lecture 15