[R-sig-Geo] Lines crossing Mercator projection map

[Giuseppe Bianco]

Thank you Roger.

Below what I get (didn't change after update.packages("rgeos"))

> library(rgeos)
rgeos: (SVN revision 357)
 GEOS runtime version: 3.3.3-CAPI-1.7.4
 Polygon checking: TRUE

>  sessionInfo()
R version 2.15.1 (2012-06-22)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] C

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets
methods   base

other attached packages:
[1] rgeos_0.2-8     rgdal_0.7-22    maptools_0.8-18 lattice_0.20-10
sp_1.0-2
[6] foreign_0.8-51

Hope this helps.

/Gis


On 12 November 2012 19:19, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> On Mon, 12 Nov 2012, Giuseppe Bianco wrote:
>
>> Dear Mike,
>>
>> Thank you for your answer.
>>
>> I get an error with the polygons intersection:
>>
>>> wrld.clip <- gIntersection(clipPoly, wrld_simpl)
>>
>>
>> Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id,
>> "rgeos_intersection")
>> :
>>
>>  TopologyException: side location conflict at -91.437497778887078
>> 17.241108073444931
>
>
> You need to provide the version of rgeos you are using - from sessionInfo(),
> and the messages printed when rgeos is loaded. For me with
>
>> library(rgeos)
>
> rgeos: (SVN revision 360)
>  GEOS runtime version: 3.3.5-CAPI-1.7.5
>  Polygon checking: TRUE
>
> and:
>
>> sessionInfo()
>
> R version 2.15.2 (2012-10-26)
> Platform: x86_64-unknown-linux-gnu (64-bit)
> ...
> other attached packages:
> [1] rgeos_0.2-9     rgdal_0.7-22    maptools_0.8-20 lattice_0.20-10
> [5] sp_1.0-2        foreign_0.8-51
>
> there is no problem.
>
>
>>
>> The problem seems to be a bad polygon in the wrld_simpl database
>>
>>> gIsValid(wrld_simpl)
>>
>>
>> [1] FALSE
>
>
> On the same system, I also see this, but it isn't the cause of your problem.
> Please also stop posting HTML, it is being stripped anyway, so post plain
> text only, as the list instructions require.
>
> Roger
>
> PS: If you need the countries, use:
>
> wrld.clip <- gIntersection(clipPoly, wrld_simpl, byid=TRUE)
> row.names(wrld.clip)
>
> shows the IDs.
>
>
>>
>> Warning message:
>>
>> In RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid") :
>>
>>  Ring Self-intersection at or near point -95.902496339999999
>> 66.946641920000005
>>
>>
>> Any way to fix it?
>>
>> BTW, I am sorry for the packages call. Here it comes:
>> library(maps)
>> library(maptools)
>> library(rgdal)
>>
>> /Gis
>>
>>
>>
>> On 12 November 2012 12:08, Michael Sumner <mdsumner at gmail.com> wrote:
>>
>>> This is happening because the maps data has data that is out of the
>>> range of [-180, 180], either from the source data itself or from the
>>> clipping somehow:
>>>
>>>  range(world.map$x, na.rm = TRUE)
>>> [1] -179.9572  190.2908
>>>
>>>
>>> You would have to carefully clean this up to make it useable, but I
>>> would just avoid it and use the clean data set in maptools (or import
>>> your own from some other source.
>>>
>>> Here you can intersect a polygon with wrld_simpl to clip it from the
>>> infinities at the poles before reprojecting:
>>>
>>> library(maptools)
>>> library(rgdal)
>>> library(rgeos)
>>>
>>> data(wrld_simpl)
>>>
>>> poly <- Polygon(cbind(c(-180, 180, 180, -180, -180), c(-80, -80, 80, 80,
>>> -80)))
>>> clipPoly <- SpatialPolygons(list(Polygons(list(poly), ID = "1")),
>>> proj4string = CRS(proj4string(wrld_simpl)))
>>>
>>> wrld.clip <- gIntersection(clipPoly, wrld_simpl)
>>>
>>> wrld.merc <- spTransform(wrld.clip, CRS("+proj=merc"))
>>>
>>> plot(wrld.merc)
>>>
>>>
>>>
>>> Also, please declare the packages you have in use with library() or
>>> require() so that code is reproducible.
>>>
>>> Cheers, Mike.
>>>
>>> On Mon, Nov 12, 2012 at 11:06 AM, Giuseppe Bianco <giubi78 at gmail.com>
>>> wrote:
>>>>
>>>> Dear List,
>>>>
>>>> I have been struggling all the day long with this issue but I didn't
>>>> find
>>>> any workaround.
>>>>
>>>> I would like to have exactly the map that come out from the code below
>>>
>>> but
>>>>
>>>> without the ugly line crossing one side to the other.
>>>>
>>>> world.map <- map(fill=T, plot=F)
>>>> IDs <- sapply(strsplit(world.map$names, ":"), function(x) x[1])
>>>> poly.map <- map2SpatialPolygons(world.map, IDs=IDs,
>>>> proj4string=CRS("+proj=longlat"))
>>>> merc.map <- spTransform(poly.map, CRS("+proj=merc"))
>>>>
>>>> # Set the map to -80 to 80
>>>> bb <- cbind(c(-180, 180), c(-80, 80))
>>>> bb.poly <- SpatialPoints(bb, proj4string = CRS("+proj=longlat"))
>>>> bb.merc <- spTransform(bb.poly, CRS("+proj=merc"))
>>>> bb.merc <- as.data.frame(bb.merc)
>>>>
>>>> # Here the ugly plot
>>>> plot(merc.map, col="gray", border="gray", xlim=bb.merc[,1],
>>>> ylim=bb.merc[,2])
>>>>
>>>> Any help is very appreciated,
>>>> /Gis
>>>>
>>>>         [[alternative HTML version deleted]]
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> R-sig-Geo at r-project.org
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>>
>>>
>>>
>>> --
>>> Michael Sumner
>>> Hobart, Australia
>>> e-mail: mdsumner at gmail.com
>>>
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> --
> Roger Bivand
> Department of Economics, NHH Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; fax +47 55 95 95 43
> e-mail: Roger.Bivand at nhh.no
>