[R-sig-Geo] Plotting GPS Coordinates on a Polar Projection Using rgdal
Michael Sumner
mdsumner at gmail.com
Mon Oct 23 22:35:53 CEST 2017
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]]
More information about the R-sig-Geo
mailing list