Plotting GPS Coordinates on a Polar Projection Using rgdal

[lmort94]

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.

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

> ## 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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[Michael Sumner-2]

On Tue, 24 Oct 2017 at 03:42 Mortimer, Linsey Anne <
[hidden email]> 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
> [hidden email]
> 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]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[lmort94]

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 <[hidden email]>
Sent: 23 October 2017 21:35
To: Mortimer, Linsey Anne
Cc: [hidden email]
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 <[hidden email]> 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
[hidden email]
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

     
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[clark richards]

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" <[hidden email]>
> To: Michael Sumner <[hidden email]>
> Cc: "[hidden email]" <[hidden email]>
> 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 <[hidden email]>
> Sent: 23 October 2017 21:35
> To: Mortimer, Linsey Anne
> Cc: [hidden email]
> 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 <
> [hidden email]> 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
> [hidden email]
> 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]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[Michael Sumner-2]

(thanks Clark! oce is excellent)

On Tue, 24 Oct 2017, 23:29 clark richards, <[hidden email]> 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" <[hidden email]>
>> To: Michael Sumner <[hidden email]>
>> Cc: "[hidden email]" <[hidden email]>
>
>
>> Subject: Re: [R-sig-Geo] Plotting GPS Coordinates on a Polar
>>         Projection Using rgdal
>>
> Message-ID:
>>         <
>> [hidden email]
>> >
>>
>> 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 <[hidden email]>
>> Sent: 23 October 2017 21:35
>> To: Mortimer, Linsey Anne
>> Cc: [hidden email]
>> 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 <
>> [hidden email]> 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
>> [hidden email]
>> 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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo