Vjylku

How is this property called for mod?

What to do with threats of blacklisting?

How vim overwrites readonly mode?

How can I prevent an oracle who can see into the past from knowing everything that has happened?

Why some papers can be published in top-tier journals while others cannot?

How much light is too much?

Is it possible to detect 100% of SQLi with a simple regex?

What are some ways of extending a description of a scenery?

"Starve to death" Vs. "Starve to the point of death"

Case protection with emphasis in biblatex

How to deal with an underperforming subordinate?

Does the ditching switch allow an A320 to float indefinitely?

How to politely refuse in-office gym instructor for steroids and protein

How to create a label containing values from different layers in QGIS

Am I correct in stating that the study of topology is purely theoretical?

What species should be used for storage of human minds?

Is it really OK to use "because of"?

Is there any advantage in specifying './' in a for loop using a glob?

Allow console draw poker game to output more hands

Buying a "Used" Router

Reading Mishnayos without understanding

Why are samba client and NFS client used differently?

Non-Cancer terminal illness that can affect young (age 10-13) girls?

How to completely remove a package in Ubuntu (like it never existed)

R: overlay a line on a raster stack/brick, and get profile of cell values along the line

R: Download a large DEM, change projection, and adjust to smaller scaleprojection differences in RReprojecting a population raster from LAEA to lonlatR: extract values from raster stack and aggregate resultsPC and Mac treating raster brick input differently in RProcessing vector to raster faster with RWhy do I get “Error: Failure during raster IO” when extracting values from my raster stack?re: Grid Issue with environmental data-frame in RRaster alignment doesn't produce aligned resultsGet colortable from three-band RGB raster (brick) with R

6

3

I've combed through R's raster functions and vignettes and can't seem to get this working.

I want to specify a line/vector across a raster stack (a DEM and possibly related variables), and get a profile of values for the cells which the line intersects. I've been able to do something similar using mask with a polygon.

EDIT: Thanks to scw, I have developed the following solution.

# I have a stack of environmental rasters in this format

new_r <- raster(ncol=615, nrow=626, xmn=-156.2, xmx=-154.8, ymn=18.89, ymx=20.30)

res(new_r) <- 0.00225

projection(new_r) <- "+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0"

values(new_r) <- outer(seq_len(nrow(new_r)), seq_len(ncol(new_r)), "+")

stackdata <- stack(new_r, sqrt(new_r))



# I designate two transect lines by long/lat

cds1 <- rbind(c(-156, 19), c(-155.5, 20.2))

cds2 <- rbind(c(-155, 20.2), c(-155, 19.2))

transects <- SpatialLines(list(Lines(list(Line(cds1)), ID = "one"), 

                               Lines(list(Line(cds2)), ID = "two")))



# plot the lines to confirm placement

plot(new_r)

plot(transects, add = TRUE)



# and return a list whose length is equal to the number of line segments,

# and each list element is a matrix with a column for each raster layer

e <- extract(stackdata, transects)

raster r

share|improve this question

edited Feb 22 '11 at 21:28

J. Win.

asked Feb 22 '11 at 8:48

J. Win.J. Win.

248414

• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  extract function prints an annoying long string of 1s while it runs. This can be hidden with invisible(capture.output(e <- extract(...))), but is there an easier way?
  
  – J. Win.
  Feb 22 '11 at 20:32
  
  
  

  
  
  
  

add a comment | 

6

3

I've combed through R's raster functions and vignettes and can't seem to get this working.

I want to specify a line/vector across a raster stack (a DEM and possibly related variables), and get a profile of values for the cells which the line intersects. I've been able to do something similar using mask with a polygon.

EDIT: Thanks to scw, I have developed the following solution.

# I have a stack of environmental rasters in this format

new_r <- raster(ncol=615, nrow=626, xmn=-156.2, xmx=-154.8, ymn=18.89, ymx=20.30)

res(new_r) <- 0.00225

projection(new_r) <- "+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0"

values(new_r) <- outer(seq_len(nrow(new_r)), seq_len(ncol(new_r)), "+")

stackdata <- stack(new_r, sqrt(new_r))



# I designate two transect lines by long/lat

cds1 <- rbind(c(-156, 19), c(-155.5, 20.2))

cds2 <- rbind(c(-155, 20.2), c(-155, 19.2))

transects <- SpatialLines(list(Lines(list(Line(cds1)), ID = "one"), 

                               Lines(list(Line(cds2)), ID = "two")))



# plot the lines to confirm placement

plot(new_r)

plot(transects, add = TRUE)



# and return a list whose length is equal to the number of line segments,

# and each list element is a matrix with a column for each raster layer

e <- extract(stackdata, transects)

raster r

share|improve this question

edited Feb 22 '11 at 21:28

J. Win.

asked Feb 22 '11 at 8:48

J. Win.J. Win.

248414

• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  extract function prints an annoying long string of 1s while it runs. This can be hidden with invisible(capture.output(e <- extract(...))), but is there an easier way?
  
  – J. Win.
  Feb 22 '11 at 20:32
  
  
  

  
  
  
  

add a comment | 

I've combed through R's raster functions and vignettes and can't seem to get this working.

I want to specify a line/vector across a raster stack (a DEM and possibly related variables), and get a profile of values for the cells which the line intersects. I've been able to do something similar using mask with a polygon.

EDIT: Thanks to scw, I have developed the following solution.

# I have a stack of environmental rasters in this format

new_r <- raster(ncol=615, nrow=626, xmn=-156.2, xmx=-154.8, ymn=18.89, ymx=20.30)

res(new_r) <- 0.00225

projection(new_r) <- "+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0"

values(new_r) <- outer(seq_len(nrow(new_r)), seq_len(ncol(new_r)), "+")

stackdata <- stack(new_r, sqrt(new_r))



# I designate two transect lines by long/lat

cds1 <- rbind(c(-156, 19), c(-155.5, 20.2))

cds2 <- rbind(c(-155, 20.2), c(-155, 19.2))

transects <- SpatialLines(list(Lines(list(Line(cds1)), ID = "one"), 

                               Lines(list(Line(cds2)), ID = "two")))



# plot the lines to confirm placement

plot(new_r)

plot(transects, add = TRUE)



# and return a list whose length is equal to the number of line segments,

# and each list element is a matrix with a column for each raster layer

e <- extract(stackdata, transects)

raster r

share|improve this question

edited Feb 22 '11 at 21:28

J. Win.

asked Feb 22 '11 at 8:48

J. Win.J. Win.

248414

# I have a stack of environmental rasters in this format

new_r <- raster(ncol=615, nrow=626, xmn=-156.2, xmx=-154.8, ymn=18.89, ymx=20.30)

res(new_r) <- 0.00225

projection(new_r) <- "+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0"

values(new_r) <- outer(seq_len(nrow(new_r)), seq_len(ncol(new_r)), "+")

stackdata <- stack(new_r, sqrt(new_r))



# I designate two transect lines by long/lat

cds1 <- rbind(c(-156, 19), c(-155.5, 20.2))

cds2 <- rbind(c(-155, 20.2), c(-155, 19.2))

transects <- SpatialLines(list(Lines(list(Line(cds1)), ID = "one"), 

                               Lines(list(Line(cds2)), ID = "two")))



# plot the lines to confirm placement

plot(new_r)

plot(transects, add = TRUE)



# and return a list whose length is equal to the number of line segments,

# and each list element is a matrix with a column for each raster layer

e <- extract(stackdata, transects)

share|improve this question

edited Feb 22 '11 at 21:28

J. Win.

asked Feb 22 '11 at 8:48

J. Win.J. Win.

248414

share|improve this question

edited Feb 22 '11 at 21:28

J. Win.

asked Feb 22 '11 at 8:48

J. Win.J. Win.

248414

asked Feb 22 '11 at 8:48

J. Win.J. Win.

248414

asked Feb 22 '11 at 8:48

J. Win.J. Win.

248414

2 Answers
2

active

oldest

votes

8

The extract should do the trick, but you may need to update to the version of raster on CRAN first. To use it, pull in the geometries you're interested in into SpatialLines objects like so:

require raster

require rgdal



r <- raster('dem.tif')

lines <- readOGR(dsn='lines.shp', layer='lines')



elevations <- extract(r, lines)

This works well for most analysis, but isn't fast enough if you're performing very large sets of data (I have an OGR/GDAL implementation I can post somewhere if it'd be useful).

share|improve this answer

answered Feb 22 '11 at 9:32

scwscw

14.5k65296

• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  +1. But approximately how large is "very large"? And how fast is "fast enough"?
  
  – whuber♦
  Feb 22 '11 at 15:09
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  Thanks, I have done it slightly different but you put me on the right track. This takes under a minute with my sample data, so is "fast enough." My related question is at gis.stackexchange.com/questions/6424/…
  
  – J. Win.
  Feb 22 '11 at 20:53
  

  
  
  
  


• 
  
  
  

  
  1
  

  
  
  
  
  
  

  
  
  @whuber: I had 7.1M lines (a distance matrix between geometries), and the raster extract function without optimization was doing about 8 lines/sec or 42k points/sec or 10 days for the dataset. The OGR/GDAL version computed the entire dataset in about 4.5hrs.
  
  – scw
  Feb 22 '11 at 22:53
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  @whuber: raster is certainly an evolving package, a coworker found an area where a small optimization to the package generated a 100x speed increase.
  
  – scw
  Feb 22 '11 at 23:19
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  My faster Python based version, for anyone needing to do many extraction operations: github.com/scw/topographic-distance/blob/master/lines.py
  
  – scw
  Mar 5 '12 at 7:34
  

  
  
  
  

 | 
show 1 more comment

3

If speed is an issue, consider using RSAGA with the profiles from lines module. http://www.saga-gis.org/saga_tool_doc/7.2.0/ta_profiles_4.html

share|improve this answer

edited 6 mins ago

answered Feb 22 '11 at 22:00

johanvdwjohanvdw

5,7601839

• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  broken link, please update.
  
  – Mouad_S
  13 hours ago
  

  
  
  
  

add a comment | 

Your Answer

StackExchange.ready(function() {
var channelOptions = {
tags: "".split(" "),
id: "79"
};
initTagRenderer("".split(" "), "".split(" "), channelOptions);

StackExchange.using("externalEditor", function() {
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled) {
StackExchange.using("snippets", function() {
createEditor();
});
}
else {
createEditor();
}
});

function createEditor() {
StackExchange.prepareEditor({
heartbeatType: 'answer',
autoActivateHeartbeat: false,
convertImagesToLinks: false,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: null,
bindNavPrevention: true,
postfix: "",
imageUploader: {
brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
allowUrls: true
},
onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
});


}
});

draft saved

draft discarded

Sign up or log in

StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});

Sign up using Google

Sign up using Facebook

Sign up using Email and Password

Post as a guest

Name

Email

Required, but never shown

StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fgis.stackexchange.com%2fquestions%2f6391%2fr-overlay-a-line-on-a-raster-stack-brick-and-get-profile-of-cell-values-along%23new-answer', 'question_page');
}
);

Post as a guest

Name

Email

Required, but never shown

8

The extract should do the trick, but you may need to update to the version of raster on CRAN first. To use it, pull in the geometries you're interested in into SpatialLines objects like so:

require raster

require rgdal



r <- raster('dem.tif')

lines <- readOGR(dsn='lines.shp', layer='lines')



elevations <- extract(r, lines)

This works well for most analysis, but isn't fast enough if you're performing very large sets of data (I have an OGR/GDAL implementation I can post somewhere if it'd be useful).

share|improve this answer

answered Feb 22 '11 at 9:32

scwscw

14.5k65296

• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  +1. But approximately how large is "very large"? And how fast is "fast enough"?
  
  – whuber♦
  Feb 22 '11 at 15:09
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  Thanks, I have done it slightly different but you put me on the right track. This takes under a minute with my sample data, so is "fast enough." My related question is at gis.stackexchange.com/questions/6424/…
  
  – J. Win.
  Feb 22 '11 at 20:53
  

  
  
  
  


• 
  
  
  

  
  1
  

  
  
  
  
  
  

  
  
  @whuber: I had 7.1M lines (a distance matrix between geometries), and the raster extract function without optimization was doing about 8 lines/sec or 42k points/sec or 10 days for the dataset. The OGR/GDAL version computed the entire dataset in about 4.5hrs.
  
  – scw
  Feb 22 '11 at 22:53
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  @whuber: raster is certainly an evolving package, a coworker found an area where a small optimization to the package generated a 100x speed increase.
  
  – scw
  Feb 22 '11 at 23:19
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  My faster Python based version, for anyone needing to do many extraction operations: github.com/scw/topographic-distance/blob/master/lines.py
  
  – scw
  Mar 5 '12 at 7:34
  

  
  
  
  

 | 
show 1 more comment

8

The extract should do the trick, but you may need to update to the version of raster on CRAN first. To use it, pull in the geometries you're interested in into SpatialLines objects like so:

require raster

require rgdal



r <- raster('dem.tif')

lines <- readOGR(dsn='lines.shp', layer='lines')



elevations <- extract(r, lines)

This works well for most analysis, but isn't fast enough if you're performing very large sets of data (I have an OGR/GDAL implementation I can post somewhere if it'd be useful).

share|improve this answer

answered Feb 22 '11 at 9:32

scwscw

14.5k65296

• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  +1. But approximately how large is "very large"? And how fast is "fast enough"?
  
  – whuber♦
  Feb 22 '11 at 15:09
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  Thanks, I have done it slightly different but you put me on the right track. This takes under a minute with my sample data, so is "fast enough." My related question is at gis.stackexchange.com/questions/6424/…
  
  – J. Win.
  Feb 22 '11 at 20:53
  

  
  
  
  


• 
  
  
  

  
  1
  

  
  
  
  
  
  

  
  
  @whuber: I had 7.1M lines (a distance matrix between geometries), and the raster extract function without optimization was doing about 8 lines/sec or 42k points/sec or 10 days for the dataset. The OGR/GDAL version computed the entire dataset in about 4.5hrs.
  
  – scw
  Feb 22 '11 at 22:53
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  @whuber: raster is certainly an evolving package, a coworker found an area where a small optimization to the package generated a 100x speed increase.
  
  – scw
  Feb 22 '11 at 23:19
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  My faster Python based version, for anyone needing to do many extraction operations: github.com/scw/topographic-distance/blob/master/lines.py
  
  – scw
  Mar 5 '12 at 7:34
  

  
  
  
  

 | 
show 1 more comment

The extract should do the trick, but you may need to update to the version of raster on CRAN first. To use it, pull in the geometries you're interested in into SpatialLines objects like so:

require raster

require rgdal



r <- raster('dem.tif')

lines <- readOGR(dsn='lines.shp', layer='lines')



elevations <- extract(r, lines)

This works well for most analysis, but isn't fast enough if you're performing very large sets of data (I have an OGR/GDAL implementation I can post somewhere if it'd be useful).

share|improve this answer

answered Feb 22 '11 at 9:32

scwscw

14.5k65296

require raster

require rgdal



r <- raster('dem.tif')

lines <- readOGR(dsn='lines.shp', layer='lines')



elevations <- extract(r, lines)

share|improve this answer

answered Feb 22 '11 at 9:32

scwscw

14.5k65296

answered Feb 22 '11 at 9:32

scwscw

14.5k65296

answered Feb 22 '11 at 9:32

scwscw

14.5k65296

3

If speed is an issue, consider using RSAGA with the profiles from lines module. http://www.saga-gis.org/saga_tool_doc/7.2.0/ta_profiles_4.html

share|improve this answer

edited 6 mins ago

answered Feb 22 '11 at 22:00

johanvdwjohanvdw

5,7601839

• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  broken link, please update.
  
  – Mouad_S
  13 hours ago
  

  
  
  
  

add a comment | 

3

If speed is an issue, consider using RSAGA with the profiles from lines module. http://www.saga-gis.org/saga_tool_doc/7.2.0/ta_profiles_4.html

share|improve this answer

edited 6 mins ago

answered Feb 22 '11 at 22:00

johanvdwjohanvdw

5,7601839

• 
  
  
  

  
  

  
  
  
  
  
  

  
  
  broken link, please update.
  
  – Mouad_S
  13 hours ago
  

  
  
  
  

add a comment | 

If speed is an issue, consider using RSAGA with the profiles from lines module. http://www.saga-gis.org/saga_tool_doc/7.2.0/ta_profiles_4.html

share|improve this answer

edited 6 mins ago

answered Feb 22 '11 at 22:00

johanvdwjohanvdw

5,7601839

share|improve this answer

edited 6 mins ago

answered Feb 22 '11 at 22:00

johanvdwjohanvdw

5,7601839

answered Feb 22 '11 at 22:00

johanvdwjohanvdw

5,7601839

answered Feb 22 '11 at 22:00

johanvdwjohanvdw

5,7601839