Lines crossing Mercator projection map

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]]

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

library(rgeos)

sessionInfo()

The problem seems to be a bad polygon in the wrld_simpl database

gIsValid(wrld_simpl)

[1] FALSE

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]]

--
Michael Sumner
Hobart, Australia
e-mail: mdsumner at gmail.com

	[[alternative HTML version deleted]]

library(rgeos)

 sessionInfo()

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]]

--
Michael Sumner
Hobart, Australia
e-mail: mdsumner at gmail.com

        [[alternative HTML version deleted]]

--
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