It's been a while since I've posted. I'm deep in the data analysis and modelling mode now. I've been trying to convert information I have on vehicles which is represented as lines in a shape file to raster data, and had a major success today. I thought I would share this since I could fine no reference to the problem elsewhere. I'm using R to do my GIS manipulations, particularly the packages rgeos, raster and ncdf.
So the problem is: I have line data in a shape file, but I want it as raster data (i.e. represented in pixels). I'm interested in a particular field (in my case the number of vehicles per hour on a road). This data is represented as a line, which has a field for number of vehicles on that road in an hour. There are functions in R which can rasterize data.
So for examples I could use:
rasterize(VehAMPKhrpoll_trans, q, 'Vkmhr', fun = "sum")
where VehAMPKhrpoll is my shape file which I've imported into R and reprojected.
But I found with this function that the sum over all the lines for the field of interest does not equal the sum for this field over all the pixels. This was a major problem for me. The model for the vehicle numbers already has errors, and I didn't want to compound those errors by incorrectly rasterizing my data.
So after a few days of hacking away at R code I managed to come up with this function which uses the extent from a NETCDF file which has the extent I require, and then rasterizes the line data from my shape file, by cropping the lines within a particular pixel, then extracting the lengths of each line in the pixel, then dividing the length of that line by the original length of the line to get the proportion of the line inside the pixel, and then multiplying the field for that line (i.e. number of vehicles) by the proportion, so that the value of the field for that line inside the pixel is proportional to how much of the line is in the pixel. And in this way, if I sum over all my pixels, I obtain the sum over all my lines for that particular field of interest.
The R function I used is:
## filename1 = file name and/or location of NETCDF file
## filename2 = file name and/or location of shape file (just the first part of the file names)
## field = The field of interest in the shape file
## nrow = number of rows in the raster file
## ncol = number of columns in the raster file
rasterize_sum <-function(filename1, filename2, field, nrow, ncol) {
### Function to obtain the sum of lines in a pixel for a particular field, but where the proportional value
### of the field that is allocated to a line inside the pixel is based on the length of the line which falls
### into the pixel.
require(ncdf)
require(rgeos)
require(raster)
#Read in the NETCDF file
ncid_old <- open.ncdf( filename1, write=FALSE )
#Read in the shape file
VehAMPKhrpoll <- readOGR(dsn = getwd(), layer = filename2)
VehAMPKhrpoll_trans <- spTransform(VehAMPKhrpoll, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
#Extract the lon and lat from the NETCDF file
lon <- ncid_old$dim$lon$vals
lat <- ncid_old$dim$lat$vals
## This part is here because the coordinates in a NETCDF file gives the centre value of each pixel
lon_min <- min(lon)-(lon[2]-lon[1])*2^(-1)
lon_max <- max(lon)+(lon[2]-lon[1])*2^(-1)
lat_min <- min(lat)-(lat[2]-lat[1])*2^(-1)
lat_max <- max(lat)+(lat[2]-lat[1])*2^(-1)
space_to_centre_lon <- ((lon_max - lon_min)*ncol^(-1))*2^(-1)
space_to_centre_lat <- ((lat_max - lat_min)*nrow^(-1))*2^(-1)
lon1 <- seq((lon_min + space_to_centre_lon),(lon_max - space_to_centre_lon), length = ncol)
lat1 <- seq((lat_min + space_to_centre_lat),(lat_max - space_to_centre_lat), length = nrow)
## Getting the lengths for each line in the original shape file
BB <- gLength(VehAMPKhrpoll_trans, byid = TRUE)
## A box to get the extent of the NETCDF file
dd <- cbind(x=c(lon_min, lon_min, lon_max, lon_max, lon_min), y=c(lat_min, lat_max, lat_max, lat_min, lat_min))
#Initialising the raster file
Vkmhr_RASTER <- raster(ncol = ncol, nrow = nrow)
extent(Vkmhr_RASTER) <- extent(dd)
# a loop to go through each pixel in the raster file
for (i in 1:length(lon1)) {
for (j in 1:length(lat1)) {
# Obtaining the extent for a pixel
xmin <- lon1[i] - (lon1[2]-lon1[1])/2 ; xmax <- lon1[i] + (lon1[2]-lon1[1])/2; ymin <- lat1[j]- (lat1[2]-lat1[1])/2; ymax <- lat1[j] + (lat1[2]-lat1[1])/2
aa <- cbind(x=c(xmin, xmin, xmax, xmax, xmin), y=c(ymin, ymax, ymax, ymin, ymin))
# cropping the SpatialLinesDataFrame by the current pixel
VEH1 <- crop(VehAMPKhrpoll_trans, extent(aa))
# Initialising a vector which will contain the proportions of each line in the pixel
CC <- NULL
# Initialising a vector which will contain the proportional values of the field of interest
VEHprop <- NULL
#Filling the raster file, from bottom row to top row (lowest lats first)
#Row 1 of the raster file will be for the maximum latitude value
if (is.null(VEH1)) {
values(Vkmhr_RASTER)[(nrow-j)*ncol + i] <- NA
} else {
AA <- gLength(VEH1, byid = TRUE)
for (k in 1:length(names(AA))) {
CC[k] <- AA[names(AA)==names(AA)[k]]/BB[names(BB)==names(AA)[k]]
##BAD BAD hard coded the field value here, need it to be input in function
VEHprop[k] <- VEH1@data[,field][k]*CC[k]
}
values(Vkmhr_RASTER)[(nrow-j)*ncol + i] <- sum(VEHprop, na.rm = TRUE)
}
}
}
Vkmhr_RASTER
}
filename1 <- "ccam_filestructure.nc" ## The NETCDF file
filename2 <- "VehAMPKhrpoll" ## The shape file
Vkmhr_RASTER1 <- rasterize_sum(filename1, filename2, "Vkmhr", 25, 25) ## Be sure to put the field name
## in inverted commas
 |
| Image of the rasterized shape file lines data for vehicle kilometre numbers on Western Cape roads |
Hopefully this helps someone out there!
No comments:
Post a Comment