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

Sat Oct 28 02:31:18 CEST 2017

(thanks Clark! oce is excellent)

On Tue, 24 Oct 2017, 23:29 clark richards, <clark.richards at gmail.com> wrote:

> 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
>> <https://maps.google.com/?q=203+Channel+Highway+%0D+Kingston+Tasmania+7050+Australia&entry=gmail&source=g>
>> Kingston Tasmania 7050 Australia
>> <https://maps.google.com/?q=203+Channel+Highway+%0D+Kingston+Tasmania+7050+Australia&entry=gmail&source=g>
>>
>>
>>
>>
>> --
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

	[[alternative HTML version deleted]]