{"word_count": 962, "character_count": 5414}
We recently uploaded on http://hal.archives-ouvertes.fr/hal-00725090 a revised version of our work, with Ewen Gallic (a.k.a. @3wen) on Visualizing spatial processes using Ripley's correction: an application to bodily-injury car accident locationIn this paper, we investigate (and extend) Ripley's circumference method to correct bias of density estimation of edges (or frontiers) of regions. The idea of the method was theoretical and di#cult to implement. We provide a simple technique - based of properties of Gaussian kernels - to compute e#efficiently weights to correct border bias on frontiers of the region of interest, with an automatic selection of an optimal radius for the method. An illustration on location of bodily-injury car accident (and hot spots) in the western part of France is discussed, where a lot of accident occur close to large cities, next to the sea.Sketches of the R code can be found in the paper, to produce maps, an to describe the impact of our boundary correction. For instance, in Finistère, the distribution of car accident is the following (with a standard kernel on the left, and with correction on the right), with 186 claims (involving bodily injury)and in Morbihan with 180 claims, observed in a specific year (2008 as far as I remember),The code is the same as the one mentioned last year, except perhaps plotting functions. First, one needs to define a color scale and associated breaksbreaks &#60;- seq( min( result $ZNA , na.rm = TRUE ) * 0.95 , max ( result$ZNA , na.rm = TRUE ) * 1.05 , length = 21)col &#60;- rev( heat . colors (20) )tofinally plot the estimationimage . plot ( result $X, result $Y, result $ZNA , xlim = range (pol[,1]) , ylim = range (pol[, 2]) , breaks = breaks , col = col ,xlab = "", ylab = "", xaxt = "n", yaxt = "n", bty = "n",zlim = range ( breaks ), horizontal = TRUE )It is possible to add a contour, the observations, and the border of the polygoncontour ( result $X, result $Y, result $ZNA , add = TRUE , col = "grey ")points (X[, 1], X[, 2], pch = 19, cex = 0.2 , col = " dodger blue")polygon (pol , lwd = 2)Now, if one wants to improve the aesthetics of the map, by adding a Google Maps base map, thefirst thing to do - after loading ggmap package - is to get the base maptheMap &#60;- get_map( location =c( left =min (pol [ ,1]) , bottom =min (pol[ ,2]) , right =max (pol [ ,1]) , top =max (pol [ ,2])), source =" google ", messaging =F, color ="bw")Of course, data need to be put in the right formatgetMelt &#60;- function ( smoothed ){res &#60;- melt ( smoothed $ZNA)res [ ,1] &#60;- smoothed $X[res [ ,1]]res [ ,2] &#60;- smoothed $Y[res [ ,2]]names (res) &#60;- list ("X","Y","ZNA")return (res )}smCont &#60;- getMelt ( result )Breaks and labels should be preparedtheLabels &#60;- round (breaks ,2)indLabels &#60;- floor (seq (1, length ( theLabels ),length .out =5)) indLabels [ length ( indLabels )] &#60;- length ( theLabels ) theLabels &#60;- as. character ( theLabels [ indLabels ])theLabels [ theLabels =="0"] &#60;- " 0.00 "Now, the map can be builtP &#60;- ggmap ( theMap )P &#60;- P + geom _ point (aes(x=X, y=Y, col=ZNA), alpha =.3 , data =smCont [!is.na( smCont $ZNA ) ,], na.rm= TRUE )It is possible to add a contourP &#60;- P + geom _ contour ( data = smCont [!is.na( smCont $ZNA) ,] ,aes(x=X, y=Y, z=ZNA ), alpha =0.5 , colour =" white ")and colors need to be updatedP &#60;- P + scale _ colour _ gradient ( name ="", low=" yellow ", high ="red", breaks = breaks [ indLabels ], limits = range ( breaks ),labels = theLabels )To remove the axis legends and labels, the theme should be updatedP &#60;- P + theme ( panel . grid . minor = element _ line ( colour =NA), panel. grid . minor = element _ line ( colour =NA), panel . background =element _ rect ( fill =NA , colour =NA), axis . text .x= element _ blank() ,axis . text .y= element _ blank () , axis . ticks .x= element _ blank() ,axis . ticks .y= element _ blank () , axis . title = element _ blank() , rect = element _ blank ())Thefinal step, in order to draw the border of the polygonpolDF &#60;- data . frame (pol)colnames ( polDF ) &#60;- list ("lon","lat")(P &#60;- P + geom _ polygon ( data =polDF , mapping =( aes(x=lon , y=lat)), colour =" black ", fill =NA))Then, we've applied that methodology to estimate the road network density in those two regions, in order to understand if high intensity means that it is a dangerous area, or if it simply because there is a lot of traffic (more traffic, more accident),We have been using the dataset obtained from the Geofabrik website which providesOpen-StreetMap data. Each observation is a section of a road, and contains a few points identified by their geographical coordinates that allow to draw lines. We have use those points to estimate a proxy of road intensity, with weight going from 10 (highways) to 1 (service roads).splitroad &#60;- function ( listroad , h = 0.0025) {pts = NULLweights &#60;- types . weights [ match ( unique ( listroad $ type ), types .weights $ type ), " weight "]for (i in 1:( length ( listroad ) - 1)) {d = diag (as. matrix ( dist ( listroad [[i]]))[, 2: nrow ( listroad[[i ]]) ]))}}return (pts )}See Ewen's blog for more details on the code, http://editerna.free.fr/blog/.... Note that Ewen did publish a poster of the paper (in French), for the http://r2013-lyon.sciencesconf.org/ conference, that will be organized in Lyon on June 27th-28th, seeAll comments are welcome...