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

Tue Oct 24 00:53:22 CEST 2017

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