[R-sig-Geo] Plotting GPS Coordinates on a Polar Projection Using rgdal

Tue Oct 24 14:29:16 CEST 2017

I would definitely defer to Michael and the other true "Geo" folks on this
list**, but if you're ok with base graphics (and a much less technical
interface) you could accomplish what you want using the `oce` package:

library(oce)
data(coastlineWorld)
circ <- list(longitude=seq(-180, 180), latitude=-rep(60, 361)) # fake
circumnavigation
## use mapPlot to make a polar stereographic plot of antarctica
mapPlot(coastlineWorld, projection="+proj=stere +lat_0=-90",
latitudelim=c(-90, -45), longitudelim = c(-180, 180), col='grey')
## now add the circumnavigation with mapPoints
mapPoints(circ)

Hope that helps!
Clark

**as you can see mapPlot is using proj/gdal interface under the hood, but
being oceanographers and only partially understanding the details of the
transforms, datums, projections, etc, we tried to write an interface that
will allow us to make pretty maps without getting too far into the details.

Message: 4
> Date: Mon, 23 Oct 2017 22:53:22 +0000
> From: "Mortimer, Linsey Anne" <l.mortimer.14 at aberdeen.ac.uk>
> To: Michael Sumner <mdsumner at gmail.com>
> Cc: "r-sig-geo at r-project.org" <r-sig-geo at r-project.org>
> Subject: Re: [R-sig-Geo] Plotting GPS Coordinates on a Polar
>         Projection Using rgdal
> Message-ID:
>         <VI1PR0101MB23173329C1E17653F09ECC8086460 at VI1PR0101MB2317.
> eurprd01.prod.exchangelabs.com>
>
> Content-Type: text/plain; charset="Windows-1252"
>
> This was a big help, thank you. The coordinates are plotted and look good
> except some of my GPS coordinates are outside of the plot - I should have
> check this beforehand. Is there a way to produce a plot using spbabel
> similar to this:?
>
>
> library(ggplot2)
> world <- map_data("world")
> worldmap <- ggplot(world, aes(x=long, y=lat, group=group)) +
> ? geom_path() +
> ? scale_y_continuous(breaks=(-2:2) * 30) +
> ? scale_x_continuous(breaks=(-4:4) * 45)
> worldmap + coord_map("ortho", orientation=c(-90, 0, 0))
>
>
> or would it be better to simply use ggplot for this?
>
>
> Many thanks,
>
> Linsey
>
>
> From: Michael Sumner <mdsumner at gmail.com>
> Sent: 23 October 2017 21:35
> To: Mortimer, Linsey Anne
> Cc: r-sig-geo at r-project.org
> Subject: Re: [R-sig-Geo] Plotting GPS Coordinates on a Polar Projection
> Using rgdal
> ?
>
>
>
> On Tue, 24 Oct 2017 at 03:42 Mortimer, Linsey Anne <
> l.mortimer.14 at aberdeen.ac.uk> wrote:
>   Hello,
>
> I am new to R and attempting to plot a map of Antarctica with the GPS
> coordinates from a circumpolar navigation of the Southern Ocean. I have
> searched the archives and rgdal-related manuals/blogs/papers to find
> variations of the code below. I have two issues:  firstly, removing the
> vertical line at -90 degrees on the plot of Antarctica. Secondly, adding
> the coordinates to this plot using points(). I have tried transforming the
> coordinates of the points to fit the projection but every time the output
> is a single  point plotted at the end of the vertical line on the
> Antarctica map, a separate plot of points in a world map projection
> (stretched) rather than polar projection, or simply an error about coercing
> an S4 class to a vector. Any advice on how to correct the code  or a point
> in the direction of a guide to this kind of mapping and spatial analysis
> would be of great help.
>
>
>
>
>
>
>
>
> You can remove the dummy seam at the south pole by decomposing to raw
> coordinates and reconstructing with spbabel:?
>
>
>
> data("wrld_simpl", package = "maptools")
> crs <- sp::proj4string(wrld_simpl)
> antarctica0 <- wrld_simpl[wrld_simpl$NAME == "Antarctica", ]
> library(spbabel)
> library(dplyr)
> ant <- spbabel::sptable(antarctica0) %>%?
> ? dplyr::filter(y_ > -90) %>%?
> ? spbabel::sp(crs = crs)
>
>
>
>
>
> pr <- "+proj=laea +lat_0=-90 +ellps=WGS84 +datum=WGS84 +no_defs
> +towgs84=0,0,0"
>
>
>
> antarctica <- sp::spTransform(ant, sp::CRS(pr))
>
>
>
>
> That makes sense when plotted in polar form, but will wrap improperly in
> native longitude/latitude form.?
>
>
> ?  An idea of what my dataset looks like (10774579 obs. of 9 variables):
>
> > head(GPS)
>
> ?..GPS_fix? ? GPS_date? GPS_time GPS_milliseconds? Latitude Longitude
> GPS_status
>
> 1? ? ? ? ? 0? 2016-12-24? 14:26:03? ? ? ? ? ? ? 878 -44.76846? 35.36859? ?
> ? ? ? 2
>
> 2? ? ? ? ? 1? 2016-12-24? 14:26:04? ? ? ? ? ? ? 309 -44.76846? 35.36859? ?
> ? ? ? 2
>
> 3? ? ? ? ? 2? 2016-12-24? 14:26:04? ? ? ? ? ? ? 878 -44.76849? 35.36866? ?
> ? ? ? 1
>
> 4? ? ? ? ? 3? 2016-12-24? 14:26:05? ? ? ? ? ? ? 309 -44.76849? 35.36866? ?
> ? ? ? 2
>
> 5? ? ? ? ? 4? 2016-12-24? 14:26:05? ? ? ? ? ? ? 878 -44.76852? 35.36873? ?
> ? ? ? 1
>
> 6? ? ? ? ? 5? 2016-12-24? 14:26:06? ? ? ? ? ? ? 305 -44.76852? 35.36873? ?
> ? ? ? 2
>
> GPS_filename? ? Distance
>
> 1? D:\\Leg1\\0 CT to Marion\\12500\\ACE-D20161224-T142556.raw 0.000000000
>
> 2? D:\\Leg1\\0 CT to Marion\\12500\\ACE-D20161224-T142556.raw 0.000000000
>
> 3? D:\\Leg1\\0 CT to Marion\\12500\\ACE-D20161224-T142556.raw? ? ? ? ? NA
>
> 4? D:\\Leg1\\0 CT to Marion\\12500\\ACE-D20161224-T142556.raw 0.003509715
>
> 5? D:\\Leg1\\0 CT to Marion\\12500\\ACE-D20161224-T142556.raw? ? ? ? ? NA
>
> 6? D:\\Leg1\\0 CT to Marion\\12500\\ACE-D20161224-T142556.raw 0.007045702
>
>
>
>
> The code I?m using:
>
> > install.packages(c("sp", "raster", ?crs?))
>
> > library(maptools)
>
> > gpclibPermit()
>
> > data(wrld_simpl)
>
> > antarctica <- wrld_simpl[wrld_simpl$NAME == "Antarctica", ]
>
> > library(rgdal)
>
> > pr <- "+proj=laea +lat_0=-90 +ellps=WGS84 +datum=WGS84 +no_defs
> +towgs84=0,0,0" #define projection
>
> > antarctica.laea <- spTransform(antarctica, CRS(pr))
>
> > plot(antarctica.laea)
>
> > ## Now I attempt to make the GPS table into a spatial object and plot
> the coordinates
>
> > attach(GPS)
>
> > coordinates(GPS) <- c("Longitude","Latitude")
>
> > points(Longitude, Latitude, pch=19, col="red", cex=0.5) #This only plots
> a singular point at the end of vertical line on the map
>
>
>
>
>
>
> Here you are mixing old attach() idiom with more modern sp idiom to assign
> the coordinates, it's probably best to not use attach at all.?
>
>
>
>
> coordinates(GPS) <- c("Longitude","Latitude")
>
>
>
> ## despite the name of the coordinates, we don't?
> ## otherwise know the coordinate system in use,?
> ## so assign is here:?
> proj4string(GPS) <- sp::CRS("+init=epsg:4326")
>
>
> ## now we can transform to the projection in use
> gps.ant <- sp::spTransform(GPS, sp::CRS(pr))
> ?
>
>
> You might find the sf package easier to use now, though there are dozens
> of tracking-specific packages with some other improvements. I find
> generally that ggplot2 provides the best overall approach to plotting, but
> it requires specific handling of the  coordinate systems for each data set,
> as done here.?
>
>
> The trip package has some methods to generate SpatialLinesDataFrame (and
> other high level classes), but to do this by removing track-information to
> get to GIS-friendly form. (It's also really painful to use and is in
> desperate need of a user-friendly update.)  There's no real stand out good
> approach for this general problem, so feel free to ask more specifics.?
> There are a huge number of disparate tracking-packages listed on CRAN so
> you might explore them for your needs:?
>
>
> https://cran.r-project.org/web/views/SpatioTemporal.html
>
> (sere under "Moving objects, trajectories")
>
>
> Cheers, Mike.?
>
>
>
>
>
>   > ## Also tried to transform the coordinates but receive an error
>
> >proj4string(GPS) <- CRS("+proj=laea +lat_0=-90 +ellps=WGS84 +datum=WGS84
> +no_defs +towgs84=0,0,0")
>
> > gps.ant <- spTransform(GPS, CRS(antarctica.laea))
>
> Error in CRS(antarctica.laea) :
>
> ? no method for coercing this S4 class to a vector
>
> In addition: Warning message:
>
> In is.na(projargs) :  is.na() applied to non-(list or vector) of type 'S4'
>
> > sessionInfo()
>
> R version 3.4.2 (2017-09-28)
>
> Platform: x86_64-w64-mingw32/x64 (64-bit)
>
> Running under: Windows >= 8 x64 (build 9200)
>
>
>
> Matrix products: default
>
>
>
> locale:
>
> [1] LC_COLLATE=English_United Kingdom.1252? LC_CTYPE=English_United
> Kingdom.1252
>
> [3] LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C
>
> [5] LC_TIME=English_United Kingdom.1252
>
>
>
> attached base packages:
>
> [1] stats? ? ?graphics? grDevices utils? ? ?datasets? methods? ?base
>
>
>
> other attached packages:
>
> [1] maptools_0.9-2 rgdal_1.2-13? ?sp_1.2-5
>
>
>
> loaded via a namespace (and not attached):
>
> [1] compiler_3.4.2? tools_3.4.2? ? ?foreign_0.8-69? grid_3.4.2? ? ?
> lattice_0.20-35
>
>
>
> Many thanks,
>
>
>
> Linsey Mortimer
>
>
>
> ? ? ? ? [[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
> --
>
>
> Dr. Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> 203 Channel Highway
> Kingston Tasmania 7050 Australia
>
>
>
>
>

	[[alternative HTML version deleted]]