[R-sig-Geo] Choropleth map

[Roger Bivand]

On Wed, 25 Jun 2014, Janka VANSCHOENWINKEL wrote:

> Thanks for your advice Roger!
>
> I will check you suggestion to use the sp package.

The usual solution to practical problems like yours is to go through step 
by step, and check after each step that the output is as expected (James 
Rooney suggested that the problem might be in the merge step, which is 
quite likely). Sometimes the classes of objects lead to dispatch to 
another method than one envisaged, so that the object returned does not 
have the expected structure.

For obvious reasons, fortify depends on having a well-formed input 
object, so checking at each step is prudent. I'd suggest that doing merge 
manually (in steps involving match(), and checking that the IDs do match 
at the correct lengths) is likely to help, as vector extension can be 
confusing - that in some situations R cycles too short vectors to match 
the required length.

Roger

>
> Next time, I'll also make sure to provide downloadable data, without html.
>
>
>
>
> 2014-06-25 14:00 GMT+02:00 Roger Bivand <Roger.Bivand at nhh.no>:
>
>> Please do *not* post in HTML - all R lists require plain text.
>>
>> You have not provided a reproducible example or access to downloadable
>> data with a script showing your problem. Nobody will try to copy & paste
>> from garbled HTML.
>>
>> You are using ggplot, which mangles spatial objects. If the choropleth map
>> can be created with plot or spplot methods from sp, you can go on to use
>> ggplot, but you will most likely find it harder to debug all the way to
>> ggplot without debugging the objects - so try spplot first.
>>
>> Once that works, you are free to go on to ggplot. Have you looked at Oscar
>> Perpigñán's book on displaying spatial data - your library should have it?
>> Hadley Wickham's approach is to make all the fish into fish soup (fortify),
>> from which it is hard to extract the fish for debugging. Debug the fish
>> first before making your soup, especially checking dimensions.
>>
>> Roger
>>
>>
>> On Tue, 24 Jun 2014, Janka VANSCHOENWINKEL wrote:
>>
>>  Dear list members,
>>>
>>>
>>> I am trying to make a basic choropleth map for all nuts3 regions in Europe
>>> ( = a map where
>>>
>>> regions are colored by the value of a certain variable ??? think for
>>> instance
>>>
>>> about a populations
>>>
>>> density map).
>>>
>>>
>>> I am struggling with this for over a week by now and I always get the same
>>> error message:
>>>
>>> Error: Aesthetics must either be length one, or the same length as
>>>
>>> the dataProblems:ValueMapDf$NUTS_ID
>>>
>>>
>>>
>>> I checked the book ???ggplot2??? from Hadley Weckham and several internet
>>> sources.  They
>>>
>>> document the different steps very properly, but still I don???t know what
>>> I
>>>
>>> am doing wrong.
>>>
>>>
>>> Could somebody help me out, please? Thanks in advance!
>>>
>>> *Below you can find my code. I hope I gave sufficient information,
>>> otherwise, just let me know!*
>>>
>>>
>>>
>>> library(maptools)
>>>
>>> library(ggplot2)
>>>
>>> library(ggmap)
>>>
>>>
>>> # read administrative boundaries (I call the SpatialPolygonsDataFrame:
>>> "nuts" since it contains all nuts3 regions of Europe) )
>>>
>>> nuts <- readShapePoly(fn="NUTS_RG_60M_2006")
>>>
>>> names(nuts)
>>>
>>> [1] "OBJECTID"   "NUTS_ID"    "STAT_LEVL_" "AREA"       "LEN"
>>> "Shape_Leng" "Shape_Area"
>>>
>>>
>>> # Then I have a dataset with over 60000 observations (several observations
>>> per nuts3 region). I take the average of the different
>>>
>>> observations per nuts3 region and save it as "Value"
>>>
>>> # The dataset is called Alldata and MEt is the variable I would like to
>>> map
>>> over the different nuts3 regions
>>>
>>>
>>> Value <- aggregate(Alldata$MEt,list(Alldata$nuts3),mean)
>>>
>>> colnames(Value) <- c("NUTS_ID","ValueMEt")
>>>
>>>
>>> head(Value)
>>>
>>>  NUTS_ID   ValueMEt
>>>
>>> 1   AT111 0.04856170
>>>
>>> 2   AT112 0.05448523
>>>
>>> 3   AT113 0.04563415
>>>
>>> 4   AT121 0.02164469
>>>
>>> 5   AT122 0.02664611
>>>
>>> 6   AT123 0.03163034
>>>
>>>
>>>
>>> # merge map and data: I call the result "ValueMapDf"
>>>
>>> ValueMapDf <- merge(nuts, Value, by.x="NUTS_ID", by.y="NUTS_ID")
>>>
>>> ValueMapDf <- ValueMapDf[order(ValueMapDf$NUTS_ID),]
>>>
>>>
>>>
>>> # then I dropped some regions of the SpatialPolygonsDataFrame, for which I
>>> don't have observations (I did this manually since I didn't
>>>
>>> know if there was a way to do it automatically)
>>>
>>> ValueMapDf<-ValueMapDf[ValueMapDf$OBJECTID > 441,]
>>>
>>> ValueMapDf<-ValueMapDf[ValueMapDf$OBJECTID != 444,]
>>>
>>> ValueMapDf<-ValueMapDf[ValueMapDf$OBJECTID != 447,]
>>>
>>> ???.
>>>
>>>
>>> ValueMapDf<-ValueMapDf[ValueMapDf$OBJECTID != 1927,]
>>>
>>>
>>>
>>> # So far so good: now I continue with plotting my map with the values!
>>> However, as you can see below, here something goes wrong.
>>>
>>> # ggplot mapping
>>>
>>> # data layer
>>>
>>> m0 <- ggplot(data=ValueMapDf)
>>>
>>> # fill with Value data
>>>
>>> m1 <- m0 + geom_polygon(aes(x=long, y=lat, group=ValueMapDf$NUTS_ID,
>>> fill="ValueMapDf$ValueMEt")) + coord_equal()
>>>
>>> # add map with only borders (to have visible borders)
>>>
>>> m2 <- m1 + geom_path(aes(x=long, y=lat, group=ValueMapDf$NUTS_ID),
>>> color='black')
>>>
>>> m2
>>>
>>> Error: Aesthetics must either be length one, or the same length as the
>>> dataProblems:ValueMapDf$NUTS_ID
>>>
>>>
>>>
>>> I also tried:
>>>
>>> qplot(long,lat,data=ValueMapDf, group=ValueMapDf$NUTS_ID,
>>> fill=ValueMapDf$ValueMEt, geom="polygon")
>>>
>>> Regions defined for each Polygons
>>>
>>> Error: Aesthetics must either be length one, or the same length as the
>>> dataProblems:ValueMapDf$ValueMEt, ValueMapDf$NUTS_ID
>>>
>>>
>>>
>>>
>>>
>>>  length(ValueMapDf$OBJECTID)
>>>>
>>>
>>> [1] 1122
>>>
>>>  length(ValueMapDf)
>>>>
>>>
>>> [1] 1122
>>>
>>>  length(ValueMapDf$NUTS_ID)
>>>>
>>>
>>> [1] 1122
>>>
>>>  length(ValueMapDf$ValueMEt)
>>>>
>>>
>>> [1] 1122
>>>
>>>
>>>
>>> *When I write:*
>>>
>>>
>>> m0 <- ggplot(data=ValueMapDf)
>>>
>>> # fill with Value data
>>>
>>> m1 <- m0 + geom_polygon(aes(x=long, y=lat, group=NUTS_ID,
>>> fill="ValueMapDf$ValueMEt")) + coord_equal()
>>>
>>> # add map with only borders (to have visible borders)
>>>
>>> m2 <- m1 + geom_path(aes(x=long, y=lat, group=NUTS_ID), color='black')
>>>
>>> m2
>>>
>>> It get: Error in eval(expr, envir, enclos) : object 'NUTS_ID' not found
>>>
>>>
>>>
>>> Thanks in advance!
>>>
>>>
>>> Janka
>>>
>>>         [[alternative HTML version deleted]]
>>>
>>>
>>>
>> --
>> Roger Bivand
>> Department of Economics, Norwegian School of Economics,
>> Helleveien 30, N-5045 Bergen, Norway.
>> voice: +47 55 95 93 55; fax +47 55 95 91 00
>> e-mail: Roger.Bivand at nhh.no
>>
>

-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
e-mail: Roger.Bivand at nhh.no