[R-sig-Geo] projection of ggmap:get_map() output
Oscar Perpiñan
oscar.perpinan at upm.es
Wed May 7 15:36:48 CEST 2014
Hi again,
I forgot to include a link with an interesting comment in the
OpenStreetMap wiki about the use of these maps:
http://wiki.openstreetmap.org/wiki/Exporting_calibrated_maps#Exporting_calibrated_maps_from_OSM
This OSM user proposes to "import the raw OSM XML data". This is
possible in R thanks to the osmar package
(http://osmar.r-forge.r-project.org/), which is able to produce sp
objects.
Best,
Oscar.
-----------------------------------------------------------------
Oscar Perpiñán Lamigueiro
Dpto. Ingeniería Eléctrica (ETSIDI-UPM)
Grupo de Sistemas Fotovoltaicos (IES-UPM)
URL: http://oscarperpinan.github.io
Twitter: @oscarperpinan
2014-05-07 12:18 GMT+02:00 Agustin Lobo <alobolistas at gmail.com>:
> Edzer,
>
> Thanks for the corrections. I had moved chunks and left loading the
> packages in the wrong place.
> It is correct now. Also, I've included a link to a file with just R
> code and the link to the vector shapefile.
>
> As far as I can tell, "levelplot() + layer()"
> or "spplot() + layer()"
> do not check correspondence of the proj4string of the operands on both sides of
> the +. Actually, both commands totally ignore any CRS information.
> This is why we've got into this mess
> and must be aware of providing consistent layers.
>
> Agus
>
> On Tue, May 6, 2014 at 6:02 PM, Edzer Pebesma
> <edzer.pebesma at uni-muenster.de> wrote:
>>
>>
>> On 05/06/2014 05:33 PM, Agustin Lobo wrote:
>>> I think I've clarified this question here:
>>> http://rpubs.com/alobo/getmapCRS
>>>
>>> In short, we can get exact positioning of both layers with
>>> * plot(r); plot(v, add=TRUE)
>>> * spplot() + layer() and
>>> * levelplot() + layer()
>>> but we must be careful of CRS details.
>>>
>>> Please send comments and corrections.
>>
>> Useful!
>>
>> These rpubs read nicely, but are very discouraging to try a rerun. Did
>> you upload the R-markdown file somewhere?
>>
>> Anyway, I cannot see how the script can run when levelplot is called
>> before lattice is loaded, and also the "+" and layer() come, I believe
>> from latticeExtra which has not been loaded up to that point (2.2).
>>
>> Does the "+" operator in
>>
>> levelplot() + layer()
>>
>> or
>>
>> spplot() + layer()
>>
>> check correspondence of the proj4string of the operands on both sides of
>> the + , can it be made to do so, or do we have remember your conclusions?
>>
>>> Agus
>>>
>>> On Tue, Apr 29, 2014 at 11:40 AM, Oscar Perpiñan <oscar.perpinan at upm.es> wrote:
>>>>>> I couldn't run your code because I am not able to read the file.
>>>>> Note dsn is the directory where you have the 4 files and layer is the
>>>>> name of the layer, which is the name of the files without extension.
>>>>
>>>> My fault, sorry.
>>>>
>>>>> There has to be a reprojection on the fly because the the two layers of
>>>>> the plot are on different coordinate systems. I think it is clear now that
>>>>> ggplot uses the Google "PseudoMercator" for the matrix while expects
>>>>> all user input on epsg4326
>>>>
>>>> I am not sure what do you mean with "reprojection on the fly" but
>>>> "layer" does not make any reprojection.
>>>> In fact, it is only a wrapper around
>>>> "eval(substitute(expression(...)))" to modify easily the panel
>>>> function of a trellis graphic.
>>>> However, in my code I am using "grid.raster" with its arguments width,
>>>> height, x and y, to resize the ggmap raster and to position it in the
>>>> graphic window created with spplot.
>>>>
>>>>> I think "acceptable" depends on the scale only. At coarser scale you
>>>>> just do not see the error
>>>>> because it is relative to the extent and the result converges to the
>>>>> correct reprojection.
>>>>
>>>> I think it depends on the extent. I have tested the code again with
>>>> the administrative boundaries of Spain and some of its regions
>>>> (http://biogeo.ucdavis.edu/data/gadm2/shp/ESP_adm.zip).
>>>> It fails for the country extent (as Edzer suggested and you have found
>>>> with your file) but it works correctly for a smaller extent (code
>>>> copied below).
>>>>
>>>>> In the code you provide, how do you know NAD27 is the appropriate datum?
>>>>> There is no prj file (one of the big shortcomings of shapefiles is
>>>>> that the presence
>>>>> of crs information in the form of a prj file is not compulsory). Note
>>>>> that you define
>>>>> a datum != WGS84 and you do not define the towgs argument.
>>>>
>>>> It is an example extracted from the p.121 of the 2008 edition of
>>>> "Applied Spatial Data Analysis with R". That chapter is now available
>>>> here: http://cran.r-project.org/web/packages/maptools/vignettes/combine_maptools.pdf
>>>>
>>>> Best,
>>>>
>>>> Oscar.
>>>>
>>>> #######################################################
>>>> library(ggmap)
>>>> library(latticeExtra)
>>>> library(sp)
>>>> library(maptools)
>>>>
>>>> ## http://biogeo.ucdavis.edu/data/gadm2/shp/ESP_adm.zip
>>>> setwd('/tmp')
>>>> spain <- readShapePoly('ESP_adm/ESP_adm1')
>>>> proj4string(spain) <- CRS("+proj=longlat +datum=WGS84")
>>>>
>>>> andalucia <- spain[spain$ID_1 == 1,]
>>>> catalonia <- spain[spain$ID_1 == 6,]
>>>>
>>>> ## myPoly <- spain
>>>> myPoly <- andalucia
>>>> ## myPoly <- catalonia
>>>>
>>>> bbPoly <- bbox(myPoly)
>>>> gmap <- get_map(c(bbPoly), maptype='watercolor', source='stamen', crop=FALSE)
>>>> bbMap <- attr(gmap, 'bb')
>>>> latCenter <- with(bbMap, ll.lat + ur.lat)/2
>>>> lonCenter <- with(bbMap, ll.lon + ur.lon)/2
>>>> height <- with(bbMap, ur.lat - ll.lat)
>>>> width <- with(bbMap, ur.lon - ll.lon)
>>>>
>>>> spplot(myPoly["ID_1"], fill='transparent') +
>>>> layer(grid.raster(gmap,
>>>> x=lonCenter, y=latCenter,
>>>> width=width, height=height,
>>>> default.units='native'),
>>>> under=TRUE)
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> R-sig-Geo at r-project.org
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>> --
>> Edzer Pebesma
>> Institute for Geoinformatics (ifgi), University of Münster
>> Heisenbergstraße 2, 48149 Münster, Germany. Phone: +49 251
>> 83 33081 http://ifgi.uni-muenster.de GPG key ID 0xAC227795
>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
More information about the R-sig-Geo
mailing list