CODEWITHSUNDEEP
rcoordinatesrastercoordinate-systems
Stefano
Stefano

Reputation: 11

Plot raster with latitude and elevation instead of latitude and longitude

I am using R and I have some issues to solve this problem: I have 2 rasters (same region, same resolution, same extent, same crs):

Raster A: (latitude, longitude, values_A)

class      : RasterLayer
dimensions : 832, 541, 450112  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : -122.2458, -117.7375, 35.0625, 41.99583  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
names      : values_A
values     : -753.4842, 0  (min, max)

Raster B: (latitude, longitude, elevation)

class      : RasterLayer
dimensions : 832, 541, 450112  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : -122.2458, -117.7375, 35.0625, 41.99583  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
names      : dem
values     : 40, 4176  (min, max)

I then need to obtain a plot of the raster A as a function of elevation and latitude (x= elevation, y= latitude, and as pixel values = values_A).

Moreover, I might have the need to aggregate pixels with the same latitude and the same elevation for example with a mean function.

do you have any suggestions on how to solve this problem with R?

Thank you!

Upvotes: 0

Views: 589

Answers (2)

Robert Hijmans
Robert Hijmans

Reputation: 47146

You can set up example data like this 

library(raster)
A <- B <- raster(nrow=83, ncol=54, ext=extent(-122.2458, -117.7375, 35.0625, 41.99583))
values(A) <- rep(1:nrow(A), ncol(A))
values(B) <- 1:ncell(B)

And here is a solution

x <- cbind(elevation=values(B), latitude=yFromRow(B, 1:nrow(B)), values(A))
r <- rasterFromXYZ(x)

Or with some existing data

B <- getData("alt", country='CHE')
A <- init(B, "y")
x <- na.omit(cbind(elevation=values(B), latitude=yFromRow(B, 1:nrow(B)), A=values(A)))

You probably want to do some rounding

x[,1] <- round(x[,1], -2)
x[,2] <- round(x[,2], 1)

r <- rasterFromXYZ(x)
plot(r, asp=NA)

Upvotes: 1

Cam McMains
Cam McMains

Reputation: 51

To start us off, I generated some example data that may be similar to what you're looking at. It will be helpful to include something like this in your questions in the future so that the answer will address any quirks you might have in your particular dataset.

df.values <- data.frame(lat = seq(44, 45, 0.01),
                        long = seq(90, 91, 0.01),
                        value_A = round(rnorm(n = 101,
                                              mean = 10,
                                              sd = 5),
                                        0)
                        )

df.elevate <- data.frame(lat = seq(44, 45, 0.01),
                         long = seq(90, 91, 0.01),
                         elevation = round(rnorm(n = 101,
                                                 mean = 200,
                                                 sd = 50),
                                           0)
                         )

The next steps will be done using functions found in the dplyr package. So install it, if you haven't already, and then load it into your environment like so:

library(dplyr)

It sounds like you want to take these two separate datasets and combine them so that you can manipulate and display them together. Since both datasets have the values (elevation and value_A) collected at the latitude and longitude level, these will be your joining variables. If you're new to this kind of data manipulation, the documentation for dplyr is very helpful: (https://dplyr.tidyverse.org/reference/join.html).

df.join <- left_join(df.values, df.elevate, by = c("lat", "long"))

As for getting a summary value (like mean) for each unique elevation and latitude pair, try a grouped summary. This is again from dplyr, but I've done it the "tidy" way, using pipes. If you're not sold on pipes, I recommend you check out this article: (https://www.datacamp.com/community/tutorials/pipe-r-tutorial)

df.smooth <- df.join %>%
    group_by(lat, elevation) %>%
    summarise(mean_A = mean(value_A, na.rm = T))

Using this, you can change what kind of summary statistic you want by changing the function within summarise. For example, you might want to take only the highest value.

df.smooth <- df.join %>%
    group_by(lat, elevation) %>%
    summarise(mean_A = max(value_A, na.rm = T))

If you have any other questions, comment on this answer.

Upvotes: 0

Related Questions