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

[Roger Bivand]

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