[GIS] R – Plotting shapefile with geographical (raster) data

coordinate systemrrastershapefile

Over the last few days I have been trying to plot a shapefile over data loaded into a data.frame within R. It should be a view of radar precipitation with an overlay of the country borders from a shapefile.

The data contain per pixel values of radar reflectivity together with the lat and lon coordinates of each pixel. Therefore the data.frame contains three colums and the end of the data is shown below:

        lat     lon      precip
 33290  55.376  6.013      128
 33291  55.375  6.029      128
 33292  55.374  6.044      128
 33293  55.374  6.059      128
 33294  55.373  6.075      128
 33295  55.372  6.090      128
 33296  55.371  6.105      128
 33297  55.370  6.121      128
 33298  55.369  6.136      128
 33299  55.368  6.151      128
 33300  55.367  6.167      128
 33301  55.366  6.182      128
 33302  55.365  6.197      128
 33303  55.364  6.213      128

For each of these lat/lon coordinates I have the pixel x and y grid-coordinate so I could add those two columns to the data.frame to get someting like:

row        x        y  lat     lon        precip  
19995      394      28 55.538  6.121      128
19996      395      28 55.537  6.137      128
19997      396      28 55.536  6.152      128
19998      397      28 55.535  6.167      128
19999      398      28 55.534  6.183      128

I have no trouble loading the shapefile using:

ned.lines = readShapeLines("d:/readShapefiles/NLD_adm/NLD_adm1.shp")
plot(ned.lines)

Plotting this shapefile together with the radar data has proven quite a challenge for me as I am unskilled with working with georeferenced data in R.

I have looked amongst others at using libraries from geoR, ggplot, raster, gstat etc, but none seems to offer the solution I am looking for, but this is most likely due to a lack of understanding.

What I expect to be the most efficient method is a plot of a georeferenced raster so that the shapefile can be plotted over this raster.

The lon and lat coordinates are plotted as regular lat/lon, but when I even try to do something as simple as:

test=raster(list(x=lon,y=lat,z=precip))

I get: "Error in .local(x, ...) : data are not on a regular grid"

The lat/lon coordinates have to be changed to a correct projection, but the exact way I should handle this is still unclear to me. I found this reference., but I just cannot seem to get the data plotted in such a way that I can overlay the shapefile. When I put the precip data simply in a matrix I can easily plot it using for example levelplot, but as this matrix lacks any geo-reference I cannot join it with the shapefile.

I hope someone is able to point me in the right direction as being stuck on something that should be very simple in my mind is rather frustrating.

Edit:

As I commented below I created a raster using the pixel x and y coordinates. This avoids the problem I had with the projection. I also converted the shapefile data to a table with a lat and a lon columns which was seperated by 99999 when a new shape object began.

The final code is shown here. It definately can be written more efficiently, but it does the job for me. Now I just need to do a bit of editing to get the image looking nice.

my.data=data.frame(cbind(x.pixel,y.pixel,map.data.mm))
r2=rasterFromXYZ(my.data)
plot(r2)
borders  = read.table("D:/readShapefiles/borders.txt", header = FALSE, skip=0, sep = ",")
vector.x=0
vector.y=0
n.shapes = which(borders$V1==99999)
for(i in 1:(length(n.shapes)-1)){  
  vector.x=0
  vector.y=0
  curr.vector.y=round(borders$V1[(n.shapes[i]+1):(n.shapes[i+1]-1)], digits=3)
  curr.vector.x=round(borders$V2[(n.shapes[i]+1):(n.shapes[i+1]-1)], digits=3)
  #lines(curr.vector.x,curr.vector.y)
  for(j in 1:length(curr.vector.x)){    
    curr.xy = which.min(abs(x.lon - curr.vector.x[j])+abs(y.lat - curr.vector.y[j])) 
    vector.x[j] = x.pixel[curr.xy]
    vector.y[j] = y.pixel[curr.xy]
  }
  lines(vector.x,vector.y)
}

Edit2:

I am now able to plot the data together, but the image still looks wrong with only a small dot at grid locations instead of a nicely drawn map. Any suggestions are welcome!

Code I am using is:

my.data=data.frame(cbind(lat,lon,map.data.mm))
r2=rasterFromXYZ(my.data)
plot(r2)
plot(ned.lines, add=TRUE)

Unfortunately I am not yet able to post a image here of the result, which shows only small dots instead of an actually filled field with radar precipitation:

Best Answer

To me there are really 2 unrelated questions here.

A) What is the projection of my shapefile and how can I merge/plot it with another layer.

B) What is the form of these raster data, and why are there in a table when they could actually be a grid?

There's a lot of possibilities here and not enough to go on. I'd like to know the answers to these:

I would find out the projection of the shapefile, and you cannot do that with maptools::readShapeLines.

1) read it with rgdal and transform:

library(rgdal)
ned.lines <- readOGR("d:/readShapefiles", "NLD_adm/NLD_adm1")

## this will tell you what the projection is
## find out and tell us
 proj4string(ned.lines)

You might find this makes it work with your grid data, but there's a lot of guessing here: 

plot(spTransform(ned.lines, CRS("+proj=longlat +datum=WGS84")))

2) Find out the actual projection and get back to us, you might need this if you cannot install rgdal. 

readLines("d:/readShapefiles/NLD_adm/NLD_adm1.prj")

(Note that this is reading text lines, not geometry).

If that fails because the file doesn't exist, then there's only more guesses we can do. It looks like you don't have a grid of the radar data, so I would go back to the source and get an actual grid, or you'll need to rasterize the data from the irregular points.

Related Question