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 |
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 > > 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 |
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 |
In reply to this post by lmort94
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 |
(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> >> >> >> >> >> -- 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 |
Free forum by Nabble | Edit this page |