4

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:

  • I have been able to get a little further. I found that there is a functionc in the raster package called rasterFromXYZ which allows me to put to data in a rasterlayer. I used: r2=rasterFromXYZ(my.data) and then plotted the data using levelplot(r2) using the rasterVis package. This indeed plots the data, but as small points on a map. I would like to the points be the extend of the rastergrid so my resolution is probably wrong. I saw there is a res=c(NA,NA) option in rasterVis, but this results in a "x cell sizes are not regular". The projection is probably still a problem. – Mark Broek Sep 2 '13 at 15:32
  • 1
    Thank you for you comments. Luckily I have been able to solve the problem by simply plotting the raster data as x,y coordinates in stead of lat/lon. This way I get a full picture instead of small dots at each grid coordinate. I converter the lat/lon coordinates of the shape file to the pixel coordinates and was now able to simple plot the shapefile over the raster. A sample of the code that I used to convert the shapedata is posted below. Here n.shapes are the individual shape objects, borders a table containing the lat lon coordinates which are seperated by a 99999 where a new object starts. – Mark Broek Sep 3 '13 at 14:19
  • I might be wrong but I think you have the same problem I described in the Q&A posted in my previous comment. Have you tried it? Once more, please post a sample of your data including lat, long, x, y and data. – Oscar Perpiñán Sep 4 '13 at 8:48
  • if you only see small dots, then I guess that your cell size is not correct. did you specify the cell size? if so did you do that in meters? don`t do that but specify it in degree. – Jens Dec 14 '13 at 22:10
1

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.

  • Thanks for your reaction. I managed to solve the problem, but it definately can be done in a better way so additional help is always appreciated. For your questions: A) I already knew the projection of the shapefile. This is wsg84 and the projection string is: "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0". My problem is probably with my interpretation of the radar data. These data are in regular lat-lon. – Mark Broek Sep 3 '13 at 15:02
  • B)The actual data are in a binary file I put in a data.frame with 3 columns. I tried to put these data in a raster using rasterfromXYX(my.data.frame) and that seems to work, but when I tried to plot these data I got white image with very small dots at the exact lat/lon locations instead of a nicely filled image. This is probably to do with the grid definition, but I am too unexperienced to solve it. Like stated above I tried the res=c() option in rasterVis, but got an error stating "x cell sizes are not regular". I assume I need to convert the regular lat/lon coordinates to another projection. – Mark Broek Sep 3 '13 at 15:05
  • try r <- rasterFromXYZ(d) where d is a data.frame or matrix with columns lon,lat,precip – mdsumner Sep 3 '13 at 15:06
  • Thanks for the comment, but that is exactly what I did. Plotting the resulting raster gives me a figure with those small dots. I am now able to plot the vector over the raster, but the image is not looking as it should. I will add a screenshot to my original post. – Mark Broek Sep 4 '13 at 7:16
  • Try adding information to your question and removing noise, what are x.pixel,y.pixel etc., these are not helpful and nothing is added. – mdsumner Sep 4 '13 at 8:27

Your Answer

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