1

I'm using GRASS but, in case, I'm also familiar with R and QGIS.

I'm trying to convert a vector area map into a raster one. My problem is that the resolution of the raster is too low (and I can't change it for other reasons) so many areas disappear.

I need each raster cell covered by the vector to be converted and to associate to the cell the coverage (i.e. if the 20% of the cell is covered by the vector map its value will be 20).

2 Answers 2

1

If I were you, I would first create a higher resolution raster(say 10 times the width and 10 times the height with same extents as your target raster) and rasterize your polygon into it.

You now have a high resolution raster that you then downsample to your desired resolution. You can use a mean/sum downsample tool or, if you are comfortable in R, make your own moving window (10 by 10, shift 10 each time) sum. The sum for each 10 by 10 patch, then assigned to a 1 by 1 pixel on your final raster will be simply the proportion (conveniently out of a 100) of the area covered by the polygon.

0

Thanks user1269942. I did like you said.

I'm answering my own question only to post the R code that I made to solve the problem. It's not fast and it could be done better, but it works.

First of all, I've increased the resolution of the target raster 20 times (in both directions) and then I've rasterized the vector layer using that resolution. "vegetazione.asc" is the exported raster map in AAIGrid format with 6 row header. "scala" is the ration between rasterized map and original resolution (in this case 20). The file "progress.txt" that I write is made only to check the progress of the code.

vegetazione <- read.table("Scrivania/Vegetazione.asc",skip=6,na.strings = "255")
scala <- 20
target <- data.frame(matrix(ncol = length(vegetazione[1,])/scala,nrow = length(vegetazione[,1])/scala))
for(i in 1:(length(vegetazione[,1])/scala)){
  for(j in 1:(length(vegetazione[1,])/scala)){
    cont <- 0
    for(k in ((i-1)*scala+1):(i*scala)){
      for(l in ((j-1)*scala+1):(j*scala)){
        if(!is.na(vegetazione[k,l])){
          cont <- cont+1
        }
      }
    }
    target[i,j] <- cont
    write.table(j,"Scrivania/progress.txt",append=TRUE,quote = FALSE,row.names = FALSE,col.names = FALSE)
  }
}
target <- target/(scala*scala)

"target" is the output to be saved.

PS: in my case the "vegetazione" map has 16540*3360 cells and it takes about 3 hours to run. Maybe it can be parallelized to improve the performance

2
  • good stuff getting it working. 3 hours...I would love to dive in and help you optimize if I weren't busy. a few things to try: 1) maybe converting to a "matrix" will be faster(note slightly different syntax to access elements) 2) access elements by "slicing" instead of looping to retrieve an array that you can 'apply'-na+sum. 3) try using python... see codereview.stackexchange.com/questions/86211/… for an example of a strided filter. Mar 1, 2017 at 20:43
  • thank you for the suggestions. When I did it I was a little bit in a hurry so I chose the faster way of coding it compromising with the performance of the code. I hope this could be still useful to someone :-)
    – Dalmo1991
    Mar 2, 2017 at 22:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.