I am using r to do some spatial analysis.

Here is my project, I am using raster files (that contain information on land use) and a polyline (which is a border between two regions). I am analyzing land use change for each pixel according to the position of that pixel regarding the border. In a nutshell I want to do a spatial regression discontinuity analysis.

Here a visual figure:

enter image description here

I have been able to collect the data of each of my (17 rasters) using that code :

r1 <- raster("~path/r1.tif")
r2 <- raster("~path/r2.tif")
r17 <- raster("~path/r17.tif")

rasStack = stack(r1,r2,...,r17)
points <- rasterToPoints(rasStack,spatial=TRUE)

write.csv(points, file = "~path/points.csv")

This provides me with a csv file that contain 20 variables (point_id, lat, long, r1_value, r2_value, ..., r17_value).

I do not want to aggregate the pixels into bigger ones. Hence I have about 22 million pixels.

Here is my question: I am trying to figure out how to add to my csv file one variable with the distance between each pixel and the line. I want to find a way to get a negative distance if the pixel is on the left, and a positive distance if the pixel is on the right.

I found one alternative way to do it on QGIS, by (1) converting the border shapefile into a chain using chainage, (2) create a distance matrix for each pixel to the nearest chain, (3) create a one side buffer around the line border, (4) clip the pixels to this one side buffer, (5) merge clipped pixels to the previous one and set them as either being on the right or on the left of the border.

However I am trying to do it in r to avoid doing things manually, and because I am sure there must be one easier solution. I found some functions like get.knn but did not get the expected results.

Here an image of the line and the whole raster:

enter image description here

Here an image of how the line intersects with the raster: enter image description here

  • How wiggly is the line across the raster? How do you define "left" and "right" of it? Do you guarantee the line is a single segment with two ends that intersect with the extent of the raster? – Spacedman Feb 13 at 22:11
  • @Spacedman, the line is wiggly yes (see edited image in my question). I am not very attached to what will be defined as being on the "left" or on the "right". I am more interested in being able to simply differentiate what is on the "left" from what is on the "right" and to get the distance between each pixel and the line. Yes the line is a single segment (with many nodes) with two ends that intersect with the extent of the raster (see image). – Marcel Campion Feb 14 at 8:23
  • Does the line go further than the red area? ie does it really go to the full bounding-box extent of the data? One approach would be to intersect the line with the bounding box which would create two polygonal regions, one on either side of the line, into which you can test your points. Alternatively given your line is monotonic (ie doesn't loop back) another way is some fairly simple algebra... – Spacedman Feb 14 at 14:22
  • @Spacedman so I found a solution to define on which side of the line points (pixels) are located. I simply created a one side buffer and then clip the pixels to that buffer. What do you mean by fairly simple algebra? I looking for a way to create a distance matrix that woule provide me with the nearest distance for each pixel to the line. – Marcel Campion Feb 14 at 21:07
  • How do you create a one-sided buffer in R? To tell if a point is L or R of a monotonic vertical set of line segments, find out which segment it is level with (by comparison with Y coord of segment vertices), then work out the X value of the line at that points' Y coord. – Spacedman Feb 14 at 22:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.