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

[prev in list] [next in list] [prev in thread] [next in thread] 

List:       r-sig-geo
Subject:    Re: [R-sig-Geo] Plotting GPS Coordinates on a Polar Projection Using rgdal
From:       Michael Sumner <mdsumner () gmail ! com>
Date:       2017-10-28 0:31:18
Message-ID: CAAcGz99DRDeZon6cok8Pj8FirDnK9jLhtZ77M1EmY6GDeuJ0Jw () mail ! gmail ! com
[Download RAW message or body]

(thanks Clark! oce is excellent)

On Tue, 24 Oct 2017, 23:29 clark richards, <clark.richards@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@aberdeen.ac.uk>
> > To: Michael Sumner <mdsumner@gmail.com>
> > Cc: "r-sig-geo@r-project.org" <r-sig-geo@r-project.org>
> 
> 
> > Subject: Re: [R-sig-Geo] Plotting GPS Coordinates on a Polar
> > Projection Using rgdal
> > 
> Message-ID:
> > <
> > VI1PR0101MB23173329C1E17653F09ECC8086460@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@gmail.com>
> > Sent: 23 October 2017 21:35
> > To: Mortimer, Linsey Anne
> > Cc: r-sig-geo@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@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@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]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[prev in list] [next in list] [prev in thread] [next in thread] 