[R-sig-Geo] making wrld_simpl longitude borders something other than (-180, 180)

Fri Sep 25 09:08:55 CEST 2015

Hi all,

Thanks for all of the help!  I just wanted to follow-up and post the code
that ended up for my purposes.  (Also, this may help my future self find
the code again in a few years.)

If you want to see the graphic, see: http://phylo.wikidot.com/biogeonick

Cheers, Nick

Code:
=================
#######################################################

# Make a map of Nick Matzke's world travels
#
# Using: http://www.latlong.net/
#
#######################################################

# Plot Great Circle paths:
# https://stat.ethz.ch/pipermail/r-sig-geo/2015-June/023032.html
library(geosphere)
library(gdata)    # for read.xls
library(maptools) # for e.g. wrld_simpl
library(rgdal)      # for e.g.
library(gpclib)   # for polygon clipping
library(fields)      # for colorbar.plot

wd = "/drives/GDrive/Matzke_PhD_docs/CV/_graphics/Matzke_map/"
setwd(wd)

xlsfn = "_Matzke_where_Ive_been_v1.xlsx"
xls = read.xls(xlsfn, header=TRUE, stringsAsFactors=FALSE)
head(xls)
tail(xls)

long_offset = -70

pdffn = "Matzke_map_v1.pdf"
pdf(file=pdffn, width=12, height=6)

data(wrld_simpl) #The world as a SpatialPolygonsDataFrame
#To avoid the lines crossing the map after reprojection we need to cut
the polygons at the new break:
w <- nowrapRecenter(wrld_simpl, offset = 180+long_offset, avoidGEOS=TRUE)
#Then proceed with the reprojection (note the proj4 string for a
mollweide projection with 150°E as the new center)
CRS_string = paste0("+proj=robin +lon_0=", long_offset)
wrld_reCentered <- spTransform(w, CRS(CRS_string))
#wrld_reCentered2 <- spTransform(wrld_reCentered, CRS("+proj=longlat
+lon_0=90"))

# Color countries by visiting
country_names = setNames(rep("white", length(wrld_simpl$NAME)), wrld_simpl$NAME)
country_colors <- setNames(rep("white", length(wrld_simpl$NAME)),
wrld_simpl$NAME)
countries_visited = c("Germany",
"United States",
"Canada",
"Mexico",
"United Kingdom",
"France",
"Zimbabwe",
"Zambia",
"United Republic of Tanzania",
"Kenya",
"Greece",
"South Africa",
"Namibia",
"Lesotho",
"New Zealand",
"French Guiana",
"Brazil",
"Australia",
"Swaziland",
"Botswana",
"Finland",
"Turkey",
"Angola",
"Madagascar",
"China",
"Spain",
"Chile")

country_colors[countries_visited] = rep("grey70",
times=length(countries_visited))

plot(wrld_reCentered, border="grey40", lwd=0.5, col=country_colors)

xls_original = xls
xls$latitude = jitter(xls$latitude)
xls$longitude = jitter(xls$longitude)

cols_for_lines = rainbow(n=nrow(xls))
for (i in 2:nrow(xls))
    {
    point1 = c(xls$longitude[i-1], xls$latitude[i-1])
    point2 = c(xls$longitude[i], xls$latitude[i])
    names(point1) = c("x","y")
    names(point2) = c("x","y")

    path = gcIntermediate(p1=point1, p2=point2, n=200,
addStartEnd=TRUE, breakAtDateLine=TRUE)

    if (is.null(dim(path)) == FALSE)
        {
        L1 = Line(path)
        Ls1 = Lines(L1, ID="a")
        SL1 = SpatialLines(list(Ls1), proj4string=CRS("+proj=longlat"))
        SL1
        #plot(SL1, add=TRUE)

        lw <- nowrapSpatialLines(SL1, offset = 180+long_offset)
        #Then proceed with the reprojection (note the proj4 string for
a mollweide projection with 150°E as the new center)
        lwrld_reCentered <- spTransform(lw, CRS(CRS_string))
        plot(lwrld_reCentered, add=TRUE, col=cols_for_lines[i], lwd=2)

        #lines(sppoints, col=cols_for_lines[i])
        #path = coordinates(path)
        #sppoints = SpatialPoints(coords=path,
proj4string=CRS("+proj=robin +lon_0=-80"))
        #points(sppoints, col=cols_for_lines[i], add=TRUE)
        }

    if (is.null(dim(path)) == TRUE)
        {
        for (j in 1:length(path))
            {
            path2 = path[[j]]
            #sppoints = SpatialPoints(coords=path2,
proj4string=CRS("+proj=robin +lon_0=-80"))
            #points(sppoints, col=cols_for_lines[i], add=TRUE)
            #lines(sppoints, col=cols_for_lines[i])
            L1 = Line(path2)
            Ls1 = Lines(L1, ID="a")
            SL1 = SpatialLines(list(Ls1), proj4string=CRS("+proj=longlat"))
            SL1
            lw <- nowrapSpatialLines(SL1, offset = 180+long_offset)
            #Then proceed with the reprojection (note the proj4 string
for a mollweide projection with 150°E as the new center)
            lwrld_reCentered <- spTransform(lw, CRS(CRS_string))
            plot(lwrld_reCentered, add=TRUE, col=cols_for_lines[i], lwd=2)
            }
        }
    } # END for (i in 2:nrow(xls))

plot(0,0, pch=".", col="white", xlim=c(0,1), ylim=c(0,1))
stripvals <- ( 1:100)
colorbar.plot(x=0.5, y=0.5, strip=y, horizontal=FALSE)

dev.off()
cmdstr = paste0("open ", pdffn)
system(cmdstr)

=================

On Mon, Sep 14, 2015 at 6:58 PM, englishchristophera at gmail.com <
englishchristophera at gmail.com> wrote:

>
>
> There is the minus sign part to enter. Though knowing what you are
> ultimately seeking would help clarify. That being said- dancing around the
> dateline is conspicuous for the difficulties that can present.Chris
> ------ Original message------From: Nick MatzkeDate: Mon, Sep 14, 2015 7:47
> AMTo: Frede Aakmann Tøgersen;Cc: r-sig-geo at stat.math.ethz.ch;Subject:Re:
> [R-sig-Geo] making wrld_simpl longitude borders something other than (-180,
> 180)
> Apologies, I got my wires crossed.What I meant was, how would I make a map
> where the left edge was atlongitude 100 (e.g. Thailand), and the right edge
> is at, say 99(again Thailand).  The Pacific Ocean, and the International
> Dateline, wouldbe in the middle of this map.I tried the obvious
> thing:plot(wrld_simpl, xlim = c(100, 99))...but that produces a flipped
> map...Cheers!NickOn Mon, Sep 14, 2015 at 3:38 PM, Frede Aakmann Tøgersen
> wrote:> Hi Nick>> To be honest I do not understand your question so I am
> guessing that you> probably only need to set the xlim to c(-100, 100). Here
> is an example:>> > data(wrld_simpl)> > plot(wrld_simpl)> > plot(wrld_simpl,
> xlim = c(-100, 100))>> Why would one wants to redefine the longitudes from
> the interval (-180,> 180) to (-100, 100)?>>> Yours sincerely / Med venlig
> hilsen>> Frede Aakmann Tøgersen> Specialist, M.Sc., Ph.D.> Plant
> Performance & Modeling>> Technology & Service Solutions> T +45 9730 5135> M
> +45 2547 6050> frtog at vestas.com> http://www.vestas.com>> Company reg.
> name: Vestas Wind Systems A/S> This e-mail is subject to our e-mail
> disclaimer statement.> Please refer to www.vestas.com/legal/notice> If
> you have received this e-mail in error please contact the sender.>>>
> -----Original Message-----> From: R-sig-Geo [mailto:
> r-sig-geo-bounces at r-project.org] On Behalf Of> Nick Matzke> Sent: 14.
> september 2015 06:51> To: r-sig-geo at stat.math.ethz.ch> Subject:
> [R-sig-Geo] making wrld_simpl longitude borders something other> than
> (-180, 180)>> Hi all,>> I would like to plot wrld_simpl, but have the left
> and right edges of the> map be at something like 100 degrees longitude,
> rather than (-180, 180).>> I tried messing with the bounding box with no
> luck.>> I can imagine something like burrowing into the object and
> subtracting 80> degrees from all the longitudes, but I suspect there must
> be a better way?>> Any help appreciated!> Nick>>         [[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>   [[alternative HTML
> version deleted]]_______________________________________________R-sig-Geo
> mailing listR-sig-Geo at r-project.orghttps://
> stat.ethz.ch/mailman/listinfo/r-sig-geo
>         [[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
>

	[[alternative HTML version deleted]]