[R-sig-Geo] Turning SpatialPointsDataFrame to Raster using rasterize

[Jon Olav Skoien]

Lionel,

You are not providing a reproducible example (y70 is not given), and you 
are not telling the versions of R and the packages you are using, so I 
am not sure if my comments here will solve your case.

1. The issue with the lines is difficult to answer without knowing more, 
but it could be a version issue. At least one of the older R-versions (I 
think it was one of the 2.14-versions) tended to add white lines to many 
gridded plots. Try with a more current version if this could be your case.

2. For the disappearing points, you might end up with more than 
observation location in some of the raster cells. rasterize will compute 
the cell value from the fun argument, and rasterToPoints will only 
return one point back. In the example below, 10 rows of meuse have been 
duplicated before the second call to rasterize and rasterToPoints. You 
can see that it returns only the original 155 points in this case. If 
meuse.grid had been coarser, also mnew would have fewer point than the 
original meuse-object. You can check how many points you have within the 
distance of a grid cell with zerodist().

loadMeuse()
r = raster(meuse.grid)
ras = rasterize(meuse, r)
mnew = rasterToPoints(ras, spatial = TRUE)
summary(meuse)

meuse = rbind(meuse, meuse[1:10,])
ras = rasterize(meuse, r)
mnew2 = rasterToPoints(ras, spatial = TRUE)
summary(meuse)
summary(mnew)
summary(mnew2)

Cheers,
Jon


On 01-Feb-13 13:32, Lionel Hertzog wrote:
> Dear List,
>
> I am having troubles trying to turn a spatialpointsdataframe into a
> raster, when I use the rasterize function from the raster packages I get
> strange white lines in certain area of the map and when turning values
> back to spatialpointsdataframe some of them disappeared. Below is a code
> example
>
> #create a list of points from coordinates of grid points in europe at a
> 10 arc-minute resolution
> test_point<-data.frame(Long=y70$Long,Lat=y70$Lat,Values=rnorm(length(y70[,1])))
> coordinates(test_point)<-~Long+Lat
> proj4string(test_point)<-CRS("+proj=longlat +datum=WGS84")
> summary(test_point)
> Object of class SpatialPointsDataFrame
> Coordinates:
>     min      max
> Long -10.41667 31.91668
> Lat   34.25000 71.25001
> Is projected: FALSE
> proj4string : [+proj=longlat +datum=WGS84]
> Number of points: 31143
> Data attributes:
>     Min.   1st Qu.    Median      Mean   3rd Qu.      Max.
> -4.278000 -0.672700 -0.003040 -0.000353  0.668600  3.915000
> #create a raster from the points extent and resolution (since 60
> arc-minutes make a degree, we have a resolution of #10/60)
> r<-raster(nrows=(bbox(test_point)[2,2]-bbox(test_point)[2,1])/0.1666667,ncols=(bbox(test_point)[1,2]-bbox(test_point)[1,1])/0.1666667,xmn=bbox(test_point)[1,1],xmx=bbox(test_point)[1,2],ymn=bbox(test_point)[2,1],ymx=bbox(test_point)[2,2],crs=CRS("+proj=longlat
> +datum=WGS84"))
> #then rasterize the points
> ras<-rasterize(x=test_point,y=r)
> #then when ploting white lines appears
> plot(ras)
> #and if we retrieve the non-NA values from this raster
> new_points<-rasterToPoints(ras,spatial=TRUE)
> summary(new_points)
> Object of class SpatialPointsDataFrame
> Coordinates:
>     min      max
> x -10.33333 31.83334
> y  34.33333 71.16667
> Is projected: FALSE
> proj4string : [+proj=longlat +datum=WGS84]
> Number of points: 29373
> Data attributes:
>     Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
> 1    7515   15420   15410   23220   31140
>
> the number of points went from 31143 to 29373, what is the problem?
> Thank you,


-- 
Jon Olav Skøien
Joint Research Centre - European Commission
Institute for Environment and Sustainability (IES)
Land Resource Management Unit

Disclaimer: Views expressed in this email are those of the individual and do not necessarily represent official views of the European Commission.