Subscribe to R-sig-geo feed
This is an archive for R-sig-geo, not a forum. It is not possible to post through Nabble - you may not start a new thread nor follow up an existing thread. If you wish to post, but are not subscribed to the list through the list homepage, subscribe first (through the list homepage, not via Nabble) and post from your subscribed email address. Until 2015-06-20, subscribers could post through Nabble, but policy has been changed as too many non-subscribers misunderstood the interface.
Updated: 2 hours 57 min ago

Re: Help with simple Map of US states with predefined regions Version 2

Thu, 09/13/2018 - 15:40
I know this is not a complete solution -- and it's a very different approach -- but it should at least show you a way to reliably get states colored by region.
(I also left out Alaska and Hawaii, since the point here is how to color the regions)

require(sp)
require(rgdal)

## US Census Bureau Tiger file -- polygons of each US State
## try this URL for download
##      https://www.census.gov/cgi-bin/geo/shapefiles/index.php?year=2017&layergroup=States+%28and+equivalent%29

## unzip to working directory ( '.' )
ustf <- readOGR('.', 'tl_2017_us_state', stringsAsFactors=FALSE)

## note, the Tiger file includes 6 additional territories
dim(ustf)
## [1] 56 14

## get rid of the extra six territories  (state.name comes with R)
cus <- subset(ustf, NAME %in% state.name)

## cheap rename
cus$state <- cus$NAME
cus$abb <- cus$STUSPS

## invent ridiculous groupings of states
cus$grp <- 'a'
cus$grp[11:20] <- 'b'
cus$grp[21:30] <- 'c'
cus$grp[31:40] <- 'd'
cus$grp[41:50] <- 'e'

## assign colors to the groups
cus$color <- 'red'
cus$color[cus$grp=='b'] <- 'green'
cus$color[cus$grp=='c'] <- 'blue'
cus$color[cus$grp=='d'] <- 'brown'
cus$color[cus$grp=='e'] <- 'cyan'

## exclude Alaska, Hawaii
cus <- subset(cus, !(state %in% c('Alaska','Hawaii')))

## get rid of extraneous variables (optional)
cus <- cus[ , c('state','REGION','abb', 'grp') ]

## plot colored by regions as defined in the Census Bureau Tiger file
plot(cus, col=cus$REGION, usePolypath=FALSE)

## color "1" is black, looks bad, do this instead
plot(cus, col=as.numeric(cus$REGION)+1, usePolypath=FALSE)
text(coordinates(cus), cus$abb, col='white', cex=0.75)

## colors specified by a color variable in the data
plot(cus, col=cus$color, usePolypath=FALSE)
text(coordinates(cus), cus$abb, col='white', cex=0.75)

(my preferred graphics device does not support Polypath, but probably most others do, so one can omit usePolypath=FALSE)

-Don

--
Don MacQueen
Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062
Lab cell 925-724-7509
 
 

On 9/13/18, 5:17 AM, "R-sig-Geo on behalf of Bill Poling" <[hidden email] on behalf of [hidden email]> wrote:

    Hi,
   
    I hope someone can help me finalize this please.
   
    I am coming close to what I need using variations from two ggplot2 tutorials.
   
    This first gives me the map of the US with AK & HI but I cannot figure out how to get my 5 regions colored
   
    #https://stackoverflow.com/questions/38021188/how-to-draw-u-s-state-map-with-hi-and-ak-with-state-abbreviations-centered-us?rq=1
   
   
    library(ggplot2)
    install.packages("ggalt")
    library(ggalt)     # coord_proj
    library(albersusa) # devtools::install_github("hrbrmstr/albersusa")
    install.packages("ggthemes")
    library(ggthemes)  # theme_map
    install.packages("rgeos")
    library(rgeos)     # centroids
    library(dplyr)
   
    # composite map with AK & HI
    usa_map <- usa_composite()
   
    # calculate the centroids for each state gCentroid(usa_map, byid=TRUE) %>%
      as.data.frame() %>%
      mutate(state=usa_map@data$iso_3166_2) -> centroids
   
    # make it usable in ggplot2
    usa_map <- fortify(usa_map)
   
    View(usa_map)
    t1 <- head(usa_map,n=5)
    knitr::kable(t1, row.names=FALSE, align=c("l", "l", "r", "r", "r"))
   
    #
   
      # |long      |lat      | group| order|  region|subregion |
      # |:---------|:--------|-----:|-----:|-------:|:---------|
      # |-87.46201 |30.38968 |     1|     1| alabama|NA        |
      # |-87.48493 |30.37249 |     1|     2| alabama|NA        |
      # |-87.52503 |30.37249 |     1|     3| alabama|NA        |
      # |-87.53076 |30.33239 |     1|     4| alabama|NA        |
      # |-87.57087 |30.32665 |     1|     5| alabama|NA        |
   
    usa_map <- fortify(usa_map)
    gg <- ggplot()
    gg <- gg + geom_map(data=usa_map, map=usa_map,
                        aes(long, lat, map_id=id),
                        color="#2b2b2b", size=0.1, fill=NA)
   
    gg <- gg + geom_text(data=centroids, aes(x, y, label=state), size=2) gg <- gg + coord_proj(us_laea_proj) gg <- gg + theme_map() gg
   
   
   
   
    #************************************************************************************************************************************************************************************/
   
    This second is an alternative (however not liking AK&HI, not coming into the map like scenario one above) but also ignoring new Mexico (because recognizing a seventh field value) and I suspect it will do the same for new York and new jersey etc.. when I add them to the list.
   
    Error in scan(file = file, what = what, sep = sep, quote = quote, dec = dec,  :  line 12 did not have 6 elements
   
    When I use newmexico (all one word) it appears white in the map like the other states not in the table statement
   
    #https://stackoverflow.com/questions/38777732/r-code-to-generating-map-of-us-states-with-specific-colors
   
    library(ggplot2)
   
    read.table(text="State.Code   region            St_Abbr   Num_Estab  colors
                          1          1   alaska       Ak        13123    #f7931e
                          3          1   arizona      AZ        18053    #f7931e
                          5          1   california   CA       143937    #f7931e
                          2          1   hawaii       HI       123456    #f7931e
                          4          1   nevada       NV       654321    #f7931e
                          6          1   oregon       OR       321456    #f7931e
                          7          1   washington   WA       456123    #f7931e
                          8          2   colorado     CO       987654    #787878
                          9          2   idaho        ID       13549     #787878
                         10          2   kansas       KS       94531     #787878
                         11          2   montana      MT       456321    #787878
                         12          2   new mexico   NM     582310            #787878 <---Not liking new mexico, saying not 6
                         13          2   oklahoma     OK       214567    #787878
                         14          2   texas        TX       675421    #787878
                         15          2   utah         UT       754321    #787878
                         16          2   wyoming      WY       543124    #787878 ",
    stringsAsFactors=FALSE, header=TRUE, comment.char="") -> df
   
    usa_map1 <- map_data("state")
    t1 <- head(usa_map1,n=5)
    knitr::kable(t1, row.names=FALSE, align=c("l", "l", "r", "r", "r"))
    View(usa_map1)
    #
    #   |long      |lat      | group| order|  region|subregion |
    #   |:---------|:--------|-----:|-----:|-------:|:---------|
    #   |-87.46201 |30.38968 |     1|     1| alabama|NA        |
    #   |-87.48493 |30.37249 |     1|     2| alabama|NA        |
    #   |-87.52503 |30.37249 |     1|     3| alabama|NA        |
    #   |-87.53076 |30.33239 |     1|     4| alabama|NA        |
    #   |-87.57087 |30.32665 |     1|     5| alabama|NA        |
   
   
   
    gg <- ggplot()
    #View(gg)
    gg <- gg + geom_map(data=usa_map1, map=usa_map1,
                        aes(long, lat, map_id=region),
                        color="#2b2b2b", size=0.15, fill=NA)
   
    gg <- gg + geom_map(data=df, map=usa_map1,
                        aes(fill=colors, map_id=region),
                        color="#2b2b2b", size=0.15)
   
   
    gg <- gg + geom_text(data=centroids, aes(x, y, label=state), size=2) gg <- gg + coord_proj(us_laea_proj) gg <- gg + theme_map() gg
   
   
    gg <- gg + scale_color_identity()
    gg <- gg + coord_map("polyconic")
    gg <- gg + ggthemes::theme_map()
    gg
   
    #c( "colorado", "idaho", "kansas", "montana", "new mexico", "oklahoma","texas", "utah", "wyoming") ) #c("alaska", "arizona", "california", "hawaii", "nevada", "oregon","washington"))
   
   
   
    William H. Poling, Ph.D., MPH
   
   
   
   
    Confidentiality Notice This message is sent from Zelis. ...{{dropped:13}}
   
    _______________________________________________
    R-sig-Geo mailing list
    [hidden email]
    https://stat.ethz.ch/mailman/listinfo/r-sig-geo
   

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

Re: [R] Help with simple Map of US states with predefined regions Version 2

Thu, 09/13/2018 - 13:31
Thank you Jeff.

Cannot seem to get this to work in the fashion I want it to appear no matter how many websites and packages I investigate.

Letting it go for the moment.

Always appreciate your advice Sir!

WHP
From: Jeff Newmiller <[hidden email]>
Sent: Thursday, September 13, 2018 10:29 AM
To: [hidden email]; Bill Poling <[hidden email]>; r-help ([hidden email]) <[hidden email]>
Subject: Re: [R] Help with simple Map of US states with predefined regions Version 2

Your data appear to be in fixed format, not space-delimited (or delimited by any other special character), so you should use read.fwf to read it in rather that read.table.

?read.fwf

In the future you should try to identify where your errors are or your data don't look right and ask focused questions about that (reproducible) problem rather than spilling your whole script into an email. That is, your error occurred in the single read.table command and the rest of it was working fine.

On September 13, 2018 5:14:33 AM PDT, Bill Poling <[hidden email]<mailto:[hidden email]>> wrote:
>Hi,
>
>I hope someone can help me finalize this please.
>
>I am coming close to what I need using variations from two ggplot2
>tutorials.
>
>This first gives me the map of the US with AK & HI but I cannot figure
>out how to get my 5 regions colored
>
>#https://stackoverflow.com/questions/38021188/how-to-draw-u-s-state-map-with-hi-and-ak-with-state-abbreviations-centered-us?rq=1<https://stackoverflow.com/questions/38021188/how-to-draw-u-s-state-map-with-hi-and-ak-with-state-abbreviations-centered-us?rq=1>
>
>
>library(ggplot2)
>install.packages("ggalt")
>library(ggalt) # coord_proj
>library(albersusa) # devtools::install_github("hrbrmstr/albersusa")
>install.packages("ggthemes")
>library(ggthemes) # theme_map
>install.packages("rgeos")
>library(rgeos) # centroids
>library(dplyr)
>
># composite map with AK & HI
>usa_map <- usa_composite()
>
># calculate the centroids for each state
>gCentroid(usa_map, byid=TRUE) %>%
> as.data.frame() %>%
> mutate(state=usa_map@data$iso_3166_2) -> centroids
>
># make it usable in ggplot2
>usa_map <- fortify(usa_map)
>
>View(usa_map)
>t1 <- head(usa_map,n=5)
>knitr::kable(t1, row.names=FALSE, align=c("l", "l", "r", "r", "r"))
>
>#
>
> # |long |lat | group| order| region|subregion |
> # |:---------|:--------|-----:|-----:|-------:|:---------|
> # |-87.46201 |30.38968 | 1| 1| alabama|NA |
> # |-87.48493 |30.37249 | 1| 2| alabama|NA |
> # |-87.52503 |30.37249 | 1| 3| alabama|NA |
> # |-87.53076 |30.33239 | 1| 4| alabama|NA |
> # |-87.57087 |30.32665 | 1| 5| alabama|NA |
>
>usa_map <- fortify(usa_map)
>gg <- ggplot()
>gg <- gg + geom_map(data=usa_map, map=usa_map,
> aes(long, lat, map_id=id),
> color="#2b2b2b", size=0.1, fill=NA)
>
>gg <- gg + geom_text(data=centroids, aes(x, y, label=state), size=2)
>gg <- gg + coord_proj(us_laea_proj)
>gg <- gg + theme_map()
>gg
>
>
>
>
>#************************************************************************************************************************************************************************************/
>
>This second is an alternative (however not liking AK&HI, not coming
>into the map like scenario one above) but also ignoring new Mexico
>(because recognizing a seventh field value) and I suspect it will do
>the same for new York and new jersey etc.. when I add them to the list.
>
>Error in scan(file = file, what = what, sep = sep, quote = quote, dec =
>dec, : line 12 did not have 6 elements
>
>When I use newmexico (all one word) it appears white in the map like
>the other states not in the table statement
>
>#https://stackoverflow.com/questions/38777732/r-code-to-generating-map-of-us-states-with-specific-colors<https://stackoverflow.com/questions/38777732/r-code-to-generating-map-of-us-states-with-specific-colors>
>
>library(ggplot2)
>
>read.table(text="State.Code region St_Abbr Num_Estab
>colors
> 1 1 alaska Ak 13123 #f7931e
> 3 1 arizona AZ 18053 #f7931e
> 5 1 california CA 143937 #f7931e
> 2 1 hawaii HI 123456 #f7931e
> 4 1 nevada NV 654321 #f7931e
> 6 1 oregon OR 321456 #f7931e
> 7 1 washington WA 456123 #f7931e
> 8 2 colorado CO 987654 #787878
> 9 2 idaho ID 13549 #787878
> 10 2 kansas KS 94531 #787878
> 11 2 montana MT 456321 #787878
>12 2 new mexico NM 582310 #787878 <---Not
>liking new mexico, saying not 6
> 13 2 oklahoma OK 214567 #787878
> 14 2 texas TX 675421 #787878
> 15 2 utah UT 754321 #787878
> 16 2 wyoming WY 543124 #787878 ",
>stringsAsFactors=FALSE, header=TRUE, comment.char="") -> df
>
>usa_map1 <- map_data("state")
>t1 <- head(usa_map1,n=5)
>knitr::kable(t1, row.names=FALSE, align=c("l", "l", "r", "r", "r"))
>View(usa_map1)
>#
># |long |lat | group| order| region|subregion |
># |:---------|:--------|-----:|-----:|-------:|:---------|
># |-87.46201 |30.38968 | 1| 1| alabama|NA |
># |-87.48493 |30.37249 | 1| 2| alabama|NA |
># |-87.52503 |30.37249 | 1| 3| alabama|NA |
># |-87.53076 |30.33239 | 1| 4| alabama|NA |
># |-87.57087 |30.32665 | 1| 5| alabama|NA |
>
>
>
>gg <- ggplot()
>#View(gg)
>gg <- gg + geom_map(data=usa_map1, map=usa_map1,
> aes(long, lat, map_id=region),
> color="#2b2b2b", size=0.15, fill=NA)
>
>gg <- gg + geom_map(data=df, map=usa_map1,
> aes(fill=colors, map_id=region),
> color="#2b2b2b", size=0.15)
>
>
>gg <- gg + geom_text(data=centroids, aes(x, y, label=state), size=2)
>gg <- gg + coord_proj(us_laea_proj)
>gg <- gg + theme_map()
>gg
>
>
>gg <- gg + scale_color_identity()
>gg <- gg + coord_map("polyconic")
>gg <- gg + ggthemes::theme_map()
>gg
>
>#c( "colorado", "idaho", "kansas", "montana", "new mexico",
>"oklahoma","texas", "utah", "wyoming") )
>#c("alaska", "arizona", "california", "hawaii", "nevada",
>"oregon","washington"))
>
>
>
>William H. Poling, Ph.D., MPH
>
>
>
>
>Confidentiality Notice This message is sent from Zelis.
>...{{dropped:13}}
>
>______________________________________________
>[hidden email]<mailto:[hidden email]> mailing list -- To UNSUBSCRIBE and more, see
>https://stat.ethz.ch/mailman/listinfo/r-help<https://stat.ethz.ch/mailman/listinfo/r-help>
>PLEASE do read the posting guide
>http://www.R-project.org/posting-guide.html<http://www.R-project.org/posting-guide.html>
>and provide commented, minimal, self-contained, reproducible code.
--
Sent from my phone. Please excuse my brevity.

Confidentiality Notice This message is sent from Zelis. This transmission may contain information which is privileged and confidential and is intended for the personal and confidential use of the named recipient only. Such information may be protected by applicable State and Federal laws from this disclosure or unauthorized use. If the reader of this message is not the intended recipient, or the employee or agent responsible for delivering the message to the intended recipient, you are hereby notified that any disclosure, review, discussion, copying, or taking any action in reliance on the contents of this transmission is strictly prohibited. If you have received this transmission in error, please contact the sender immediately. Zelis, 2018.

        [[alternative HTML version deleted]]

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

Help with simple Map of US states with predefined regions Version 2

Thu, 09/13/2018 - 07:17
Hi,

I hope someone can help me finalize this please.

I am coming close to what I need using variations from two ggplot2 tutorials.

This first gives me the map of the US with AK & HI but I cannot figure out how to get my 5 regions colored

#https://stackoverflow.com/questions/38021188/how-to-draw-u-s-state-map-with-hi-and-ak-with-state-abbreviations-centered-us?rq=1


library(ggplot2)
install.packages("ggalt")
library(ggalt)     # coord_proj
library(albersusa) # devtools::install_github("hrbrmstr/albersusa")
install.packages("ggthemes")
library(ggthemes)  # theme_map
install.packages("rgeos")
library(rgeos)     # centroids
library(dplyr)

# composite map with AK & HI
usa_map <- usa_composite()

# calculate the centroids for each state gCentroid(usa_map, byid=TRUE) %>%
  as.data.frame() %>%
  mutate(state=usa_map@data$iso_3166_2) -> centroids

# make it usable in ggplot2
usa_map <- fortify(usa_map)

View(usa_map)
t1 <- head(usa_map,n=5)
knitr::kable(t1, row.names=FALSE, align=c("l", "l", "r", "r", "r"))

#

  # |long      |lat      | group| order|  region|subregion |
  # |:---------|:--------|-----:|-----:|-------:|:---------|
  # |-87.46201 |30.38968 |     1|     1| alabama|NA        |
  # |-87.48493 |30.37249 |     1|     2| alabama|NA        |
  # |-87.52503 |30.37249 |     1|     3| alabama|NA        |
  # |-87.53076 |30.33239 |     1|     4| alabama|NA        |
  # |-87.57087 |30.32665 |     1|     5| alabama|NA        |

usa_map <- fortify(usa_map)
gg <- ggplot()
gg <- gg + geom_map(data=usa_map, map=usa_map,
                    aes(long, lat, map_id=id),
                    color="#2b2b2b", size=0.1, fill=NA)

gg <- gg + geom_text(data=centroids, aes(x, y, label=state), size=2) gg <- gg + coord_proj(us_laea_proj) gg <- gg + theme_map() gg




#************************************************************************************************************************************************************************************/

This second is an alternative (however not liking AK&HI, not coming into the map like scenario one above) but also ignoring new Mexico (because recognizing a seventh field value) and I suspect it will do the same for new York and new jersey etc.. when I add them to the list.

Error in scan(file = file, what = what, sep = sep, quote = quote, dec = dec,  :  line 12 did not have 6 elements

When I use newmexico (all one word) it appears white in the map like the other states not in the table statement

#https://stackoverflow.com/questions/38777732/r-code-to-generating-map-of-us-states-with-specific-colors

library(ggplot2)

read.table(text="State.Code   region            St_Abbr   Num_Estab  colors
                      1          1   alaska       Ak        13123    #f7931e
                      3          1   arizona      AZ        18053    #f7931e
                      5          1   california   CA       143937    #f7931e
                      2          1   hawaii       HI       123456    #f7931e
                      4          1   nevada       NV       654321    #f7931e
                      6          1   oregon       OR       321456    #f7931e
                      7          1   washington   WA       456123    #f7931e
                      8          2   colorado     CO       987654    #787878
                      9          2   idaho        ID       13549     #787878
                     10          2   kansas       KS       94531     #787878
                     11          2   montana      MT       456321    #787878
                     12          2   new mexico   NM     582310            #787878 <---Not liking new mexico, saying not 6
                     13          2   oklahoma     OK       214567    #787878
                     14          2   texas        TX       675421    #787878
                     15          2   utah         UT       754321    #787878
                     16          2   wyoming      WY       543124    #787878 ",
stringsAsFactors=FALSE, header=TRUE, comment.char="") -> df

usa_map1 <- map_data("state")
t1 <- head(usa_map1,n=5)
knitr::kable(t1, row.names=FALSE, align=c("l", "l", "r", "r", "r"))
View(usa_map1)
#
#   |long      |lat      | group| order|  region|subregion |
#   |:---------|:--------|-----:|-----:|-------:|:---------|
#   |-87.46201 |30.38968 |     1|     1| alabama|NA        |
#   |-87.48493 |30.37249 |     1|     2| alabama|NA        |
#   |-87.52503 |30.37249 |     1|     3| alabama|NA        |
#   |-87.53076 |30.33239 |     1|     4| alabama|NA        |
#   |-87.57087 |30.32665 |     1|     5| alabama|NA        |



gg <- ggplot()
#View(gg)
gg <- gg + geom_map(data=usa_map1, map=usa_map1,
                    aes(long, lat, map_id=region),
                    color="#2b2b2b", size=0.15, fill=NA)

gg <- gg + geom_map(data=df, map=usa_map1,
                    aes(fill=colors, map_id=region),
                    color="#2b2b2b", size=0.15)


gg <- gg + geom_text(data=centroids, aes(x, y, label=state), size=2) gg <- gg + coord_proj(us_laea_proj) gg <- gg + theme_map() gg


gg <- gg + scale_color_identity()
gg <- gg + coord_map("polyconic")
gg <- gg + ggthemes::theme_map()
gg

#c( "colorado", "idaho", "kansas", "montana", "new mexico", "oklahoma","texas", "utah", "wyoming") ) #c("alaska", "arizona", "california", "hawaii", "nevada", "oregon","washington"))



William H. Poling, Ph.D., MPH




Confidentiality Notice This message is sent from Zelis. ...{{dropped:13}}

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

Re: help: Problem getting centroids inside lists

Wed, 09/12/2018 - 16:09
I understand, but I said "it doesn't work" because it gave an undesired
solution, but I'll try to be more explicit next time.
And I forgot to say that: s  <- 1: length(ct)

And thank you for your advice, but the answer of Vijay Lulla resolved my
problem.

Regards,
Ariel

2018-09-12 17:43 GMT-03:00 MacQueen, Don <[hidden email]>:

> On any R mailing list, whenever you say "it doesn't work", you should
> always copy exactly what R command you gave, and any error messages.
>
> First notice:
>
> > ct$a$two
> SpatialPoints:
>      lon lat
> [1,] -18  -2
> [2,] -16  50
> Coordinate Reference System (CRS) arguments: NA
>
> Then compare these:
>
> > coordinates(ct$a$two)
>      lon lat
> [1,] -18  -2
> [2,] -16  50
>
> > coordinates(slot(ct$a$two,'coords'))
>      lon lat
> [1,] -18  -2
> [2,] -16  50
>
> The results are the same. Why are you using the slot() function? There is
> no need, and it makes it more difficult to understand.
>
>
> You have assumed the mean() function will give you what you want when
> there are two points. It doesn't. Try it and see:
>
> > mean( coordinates(ct$a$two))
> [1] 3.5
>
> > mean(coordinates(slot(ct$a$two,'coords')))
>  [1] 3.5
>
> Replace the mean() function with a function that will give you the
> "centroid" of two points, however you want to define that centroid. Perhaps
> the means of the two columns? R has a function for that.
>
> -----------------
> In this bit:
>
> ctply <- lapply(X = s, FUN = function(x) sapply(X = ct[[x]],
>                                                     FUN = function(y)
>
> Using X = s is wrong because you haven't defined or created s anywhere. X
> should be ct.
>
> Also, X = ct[[x]] is wrong. When the first FUN is executed, it will be
> supplied automatically with the whole of ct$a, then ct$b. Try
>
> ctply <- lapply(X = s, FUN = function(x) sapply(X = x,
>                                                     FUN = function(y)
>
> The whole thing is easier to understand and test if you define the
> functions outside the lapply and sapply calls.
>
> -Don
> --
> Don MacQueen
> Lawrence Livermore National Laboratory
> 7000 East Ave., L-627
> Livermore, CA 94550
> 925-423-1062
> Lab cell 925-724-7509
>
>
>
> On 9/12/18, 11:42 AM, "R-sig-Geo on behalf of Ariel Fuentesdi" <
> [hidden email] on behalf of [hidden email]>
> wrote:
>
>     I went to a more general approach, that is:
>
>     ctply <- lapply(X = s, FUN = function(x) sapply(X = ct[[x]],
>               FUN = function(y) if(length(y)>1) geosphere::centroid(slot(y,
>     "coords"))
>               else sp::coordinates(slot(y, "coords"))))
>
>     But now I want to add the case when they are only two elements. The
> dataset
>     will be:
>
>     ct <- list(a = list(one = data.frame(lon = c(-180, -160, -60), lat =
> c(-20,
>     5, 0)),
>                         two = data.frame(lon = c(-18, -16), lat = c(-2,
> 50))),
>                b = list(one = data.frame(lon = c(-9, -8, -3), lat = c(-1,
> 25,
>     5)),
>                         two = data.frame(lon = c(-90), lat = c(-1))))
>     coordinates(ct$a$one) <- ~lon+lat
>     coordinates(ct$a$two) <- ~lon+lat
>     coordinates(ct$b$one) <- ~lon+lat
>     coordinates(ct$b$two) <- ~lon+lat
>
>     And I did the following but it doesn't work:
>
>     ctply <- lapply(X = s, FUN = function(x) sapply(X = ct[[x]],
>                                                     FUN = function(y)
>     if(length(y)>2) geosphere::centroid(slot(y, "coords"))
>                                                     else if (length(y) ==
> 1)
>     sp::coordinates(slot(y, "coords"))
>                                                     else
>     mean(sp::coordinates(slot(y, "coords")))))
>
>     I need the result of each element of the list will be a matrix of two
> rows
>     per column (named: one, two). How do I fix it?
>
>     Regards,
>     Ariel
>
>     2018-09-10 17:29 GMT-03:00 MacQueen, Don <[hidden email]>:
>
>     > If all of your data frames had enough points then this should work:
>     >
>     > myfun1 <- function(le) {
>     >   list(one=geosphere::centroid(coordinates(le$one)),
>     >        two=geosphere::centroid(coordinates(le$two))
>     >        )
>     > }
>     >
>     > lapply(ct, myfun1)
>     >
>     > To handle the one point case, try the following, which I think works
> with
>     > your example data, but is untested if there just two points.
>     > It also assumes the data frames are named 'one' and 'two', and will
> ignore
>     > any others. To handle other names, I think you
>     > could replace  $one  with  [[1]]  and  $two  with  [[2]]  .
>     >
>     > myfun <- function(le) {
>     >   list(one=if (length(le$one)>1) geosphere::centroid(
> coordinates(le$one))
>     > else coodinates(le$one),
>     >        two=if (length(le$two)>1) geosphere::centroid(
> coordinates(le$two))
>     > else coordinates(le$two)
>     >        )
>     > }
>     >
>     > lapply(ct, myfun)
>     >
>     >
>     > To see a little bit of what's going on, try
>     >
>     >   myfun(ct[[1]])
>     >   myfun1(ct[[1]])
>     >   myfun1(ct[[2]])
>     >   myfun(ct[[2]])
>     >
>     >
>     >
>     > I would also verify that the length method for objects of class
>     > SpatialPoints does return the number of points, as it appears to:
>     >
>     > > class(ct$a$two)
>     > [1] "SpatialPoints"
>     > attr(,"package")
>     > [1] "sp"
>     >
>     > > length(ct$a$one)
>     > [1] 3
>     > > length(ct$a$two)
>     > [1] 3
>     > > length(ct$a$two)
>     > [1] 3
>     > > length(ct$b$two)
>     > [1] 1
>     >
>     > --
>     > Don MacQueen
>     > Lawrence Livermore National Laboratory
>     > 7000 East Ave., L-627
>     > Livermore, CA 94550
>     > 925-423-1062
>     > Lab cell 925-724-7509
>     >
>     >
>     >
>     > On 9/10/18, 10:56 AM, "R-sig-Geo on behalf of Ariel Fuentesdi" <
>     > [hidden email] on behalf of
> [hidden email]>
>     > wrote:
>     >
>     >     Hi everyone,
>     >
>     >     I have a list of coordinates called "ct" and I want to extract
> the
>     >     centroids of each sublist, but it only works when it has only 3
> or more
>     >     points. In ct$b$two only has one point; in that case, I would
> rescue
>     > that
>     >     point as If it were a centroid.
>     >
>     >     This is what I did:
>     >
>     >     library(dplyr)
>     >     library(geosphere)
>     >
>     >     ct <- list(a = list(one = data.frame(lon = c(-180, -160, -60),
> lat =
>     > c(-20,
>     >     5, 0)),
>     >                         two = data.frame(lon = c(-18, -16, -6), lat =
>     > c(-2, 50,
>     >     10))),
>     >                b = list(one = data.frame(lon = c(-9, -8, -3), lat =
> c(-1,
>     > 25,
>     >     5)),
>     >                         two = data.frame(lon = c(-90), lat = c(-1))))
>     >
>     >     coordinates(ct$a$one) <- ~lon+lat
>     >     coordinates(ct$a$two) <- ~lon+lat
>     >     coordinates(ct$b$one) <- ~lon+lat
>     >     coordinates(ct$b$two) <- ~lon+lat
>     >
>     >     s <- 1:length(ct)
>     >     ctply <- list()
>     >     for (i in s){
>     >      for (j in 1:length(ct[[i]])) {
>     >       ctply[[i]][j] <- ifelse(test = lengths(ct[[i]][1]) > 2, yes =
>     > sapply(X =
>     >     ct[[i]][j],
>     >         FUN = function(y) geosphere::centroid(slot(y, "coords"))),
> no =
>     >     ct[[i]][j] )
>     >      }
>     >     }
>     >
>     >     Thanks in advance
>     >
>     >     Regards,
>     >     Ariel
>     >
>     >         [[alternative HTML version deleted]]
>     >
>     >     _______________________________________________
>     >     R-sig-Geo mailing list
>     >     [hidden email]
>     >     https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>     >
>     >
>     >
>
>         [[alternative HTML version deleted]]
>
>     _______________________________________________
>     R-sig-Geo mailing list
>     [hidden email]
>     https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
>
        [[alternative HTML version deleted]]

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

Re: help: Problem getting centroids inside lists

Wed, 09/12/2018 - 15:43
On any R mailing list, whenever you say "it doesn't work", you should always copy exactly what R command you gave, and any error messages.

First notice:

> ct$a$two
SpatialPoints:
     lon lat
[1,] -18  -2
[2,] -16  50
Coordinate Reference System (CRS) arguments: NA

Then compare these:

> coordinates(ct$a$two)
     lon lat
[1,] -18  -2
[2,] -16  50

> coordinates(slot(ct$a$two,'coords'))
     lon lat
[1,] -18  -2
[2,] -16  50

The results are the same. Why are you using the slot() function? There is no need, and it makes it more difficult to understand.


You have assumed the mean() function will give you what you want when there are two points. It doesn't. Try it and see:

> mean( coordinates(ct$a$two))
[1] 3.5

> mean(coordinates(slot(ct$a$two,'coords')))
 [1] 3.5

Replace the mean() function with a function that will give you the "centroid" of two points, however you want to define that centroid. Perhaps the means of the two columns? R has a function for that.

-----------------
In this bit:

ctply <- lapply(X = s, FUN = function(x) sapply(X = ct[[x]],
                                                    FUN = function(y)

Using X = s is wrong because you haven't defined or created s anywhere. X should be ct.

Also, X = ct[[x]] is wrong. When the first FUN is executed, it will be supplied automatically with the whole of ct$a, then ct$b. Try

ctply <- lapply(X = s, FUN = function(x) sapply(X = x,
                                                    FUN = function(y)

The whole thing is easier to understand and test if you define the functions outside the lapply and sapply calls.

-Don
--
Don MacQueen
Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062
Lab cell 925-724-7509
 
 

On 9/12/18, 11:42 AM, "R-sig-Geo on behalf of Ariel Fuentesdi" <[hidden email] on behalf of [hidden email]> wrote:

    I went to a more general approach, that is:
   
    ctply <- lapply(X = s, FUN = function(x) sapply(X = ct[[x]],
              FUN = function(y) if(length(y)>1) geosphere::centroid(slot(y,
    "coords"))
              else sp::coordinates(slot(y, "coords"))))
   
    But now I want to add the case when they are only two elements. The dataset
    will be:
   
    ct <- list(a = list(one = data.frame(lon = c(-180, -160, -60), lat = c(-20,
    5, 0)),
                        two = data.frame(lon = c(-18, -16), lat = c(-2, 50))),
               b = list(one = data.frame(lon = c(-9, -8, -3), lat = c(-1, 25,
    5)),
                        two = data.frame(lon = c(-90), lat = c(-1))))
    coordinates(ct$a$one) <- ~lon+lat
    coordinates(ct$a$two) <- ~lon+lat
    coordinates(ct$b$one) <- ~lon+lat
    coordinates(ct$b$two) <- ~lon+lat
   
    And I did the following but it doesn't work:
   
    ctply <- lapply(X = s, FUN = function(x) sapply(X = ct[[x]],
                                                    FUN = function(y)
    if(length(y)>2) geosphere::centroid(slot(y, "coords"))
                                                    else if (length(y) == 1)
    sp::coordinates(slot(y, "coords"))
                                                    else
    mean(sp::coordinates(slot(y, "coords")))))
   
    I need the result of each element of the list will be a matrix of two rows
    per column (named: one, two). How do I fix it?
   
    Regards,
    Ariel
   
    2018-09-10 17:29 GMT-03:00 MacQueen, Don <[hidden email]>:
   
    > If all of your data frames had enough points then this should work:
    >
    > myfun1 <- function(le) {
    >   list(one=geosphere::centroid(coordinates(le$one)),
    >        two=geosphere::centroid(coordinates(le$two))
    >        )
    > }
    >
    > lapply(ct, myfun1)
    >
    > To handle the one point case, try the following, which I think works with
    > your example data, but is untested if there just two points.
    > It also assumes the data frames are named 'one' and 'two', and will ignore
    > any others. To handle other names, I think you
    > could replace  $one  with  [[1]]  and  $two  with  [[2]]  .
    >
    > myfun <- function(le) {
    >   list(one=if (length(le$one)>1) geosphere::centroid(coordinates(le$one))
    > else coodinates(le$one),
    >        two=if (length(le$two)>1) geosphere::centroid(coordinates(le$two))
    > else coordinates(le$two)
    >        )
    > }
    >
    > lapply(ct, myfun)
    >
    >
    > To see a little bit of what's going on, try
    >
    >   myfun(ct[[1]])
    >   myfun1(ct[[1]])
    >   myfun1(ct[[2]])
    >   myfun(ct[[2]])
    >
    >
    >
    > I would also verify that the length method for objects of class
    > SpatialPoints does return the number of points, as it appears to:
    >
    > > class(ct$a$two)
    > [1] "SpatialPoints"
    > attr(,"package")
    > [1] "sp"
    >
    > > length(ct$a$one)
    > [1] 3
    > > length(ct$a$two)
    > [1] 3
    > > length(ct$a$two)
    > [1] 3
    > > length(ct$b$two)
    > [1] 1
    >
    > --
    > Don MacQueen
    > Lawrence Livermore National Laboratory
    > 7000 East Ave., L-627
    > Livermore, CA 94550
    > 925-423-1062
    > Lab cell 925-724-7509
    >
    >
    >
    > On 9/10/18, 10:56 AM, "R-sig-Geo on behalf of Ariel Fuentesdi" <
    > [hidden email] on behalf of [hidden email]>
    > wrote:
    >
    >     Hi everyone,
    >
    >     I have a list of coordinates called "ct" and I want to extract the
    >     centroids of each sublist, but it only works when it has only 3 or more
    >     points. In ct$b$two only has one point; in that case, I would rescue
    > that
    >     point as If it were a centroid.
    >
    >     This is what I did:
    >
    >     library(dplyr)
    >     library(geosphere)
    >
    >     ct <- list(a = list(one = data.frame(lon = c(-180, -160, -60), lat =
    > c(-20,
    >     5, 0)),
    >                         two = data.frame(lon = c(-18, -16, -6), lat =
    > c(-2, 50,
    >     10))),
    >                b = list(one = data.frame(lon = c(-9, -8, -3), lat = c(-1,
    > 25,
    >     5)),
    >                         two = data.frame(lon = c(-90), lat = c(-1))))
    >
    >     coordinates(ct$a$one) <- ~lon+lat
    >     coordinates(ct$a$two) <- ~lon+lat
    >     coordinates(ct$b$one) <- ~lon+lat
    >     coordinates(ct$b$two) <- ~lon+lat
    >
    >     s <- 1:length(ct)
    >     ctply <- list()
    >     for (i in s){
    >      for (j in 1:length(ct[[i]])) {
    >       ctply[[i]][j] <- ifelse(test = lengths(ct[[i]][1]) > 2, yes =
    > sapply(X =
    >     ct[[i]][j],
    >         FUN = function(y) geosphere::centroid(slot(y, "coords"))), no =
    >     ct[[i]][j] )
    >      }
    >     }
    >
    >     Thanks in advance
    >
    >     Regards,
    >     Ariel
    >
    >         [[alternative HTML version deleted]]
    >
    >     _______________________________________________
    >     R-sig-Geo mailing list
    >     [hidden email]
    >     https://stat.ethz.ch/mailman/listinfo/r-sig-geo
    >
    >
    >
   
    [[alternative HTML version deleted]]
   
    _______________________________________________
    R-sig-Geo mailing list
    [hidden email]
    https://stat.ethz.ch/mailman/listinfo/r-sig-geo
   

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

Re: help: Problem getting centroids inside lists

Wed, 09/12/2018 - 15:06
Maybe the following function is what you're looking for?

getcentroids <- function(x1) {
  getcentroid <- function(x) {
    coords <- slot(x, "coords")
    numrows <- dim(coords)[1]
    ret <- if(numrows == 2) {
      matrix(apply(coords,2,mean), nrow=1)
    } else {
      if(numrows == 1) {
        coords
      } else {
        geosphere::centroid(coords)
      }
    }
    ret
  }

  r <- lapply(x1, function(x) as.data.frame(getcentroid(x)))

  ret <- matrix(unlist(r), ncol=2, byrow=TRUE)
  rownames(ret) <- names(r); colnames(ret) <- c("lon", "lat")

  ret
}

Now call it as lapply(ct, getcentroids)

HTH,
Vijay.

On Wed, Sep 12, 2018 at 2:42 PM Ariel Fuentesdi <[hidden email]>
wrote:

> I went to a more general approach, that is:
>
> ctply <- lapply(X = s, FUN = function(x) sapply(X = ct[[x]],
>           FUN = function(y) if(length(y)>1) geosphere::centroid(slot(y,
> "coords"))
>           else sp::coordinates(slot(y, "coords"))))
>
> But now I want to add the case when they are only two elements. The dataset
> will be:
>
> ct <- list(a = list(one = data.frame(lon = c(-180, -160, -60), lat = c(-20,
> 5, 0)),
>                     two = data.frame(lon = c(-18, -16), lat = c(-2, 50))),
>            b = list(one = data.frame(lon = c(-9, -8, -3), lat = c(-1, 25,
> 5)),
>                     two = data.frame(lon = c(-90), lat = c(-1))))
> coordinates(ct$a$one) <- ~lon+lat
> coordinates(ct$a$two) <- ~lon+lat
> coordinates(ct$b$one) <- ~lon+lat
> coordinates(ct$b$two) <- ~lon+lat
>
> And I did the following but it doesn't work:
>
> ctply <- lapply(X = s, FUN = function(x) sapply(X = ct[[x]],
>                                                 FUN = function(y)
> if(length(y)>2) geosphere::centroid(slot(y, "coords"))
>                                                 else if (length(y) == 1)
> sp::coordinates(slot(y, "coords"))
>                                                 else
> mean(sp::coordinates(slot(y, "coords")))))
>
> I need the result of each element of the list will be a matrix of two rows
> per column (named: one, two). How do I fix it?
>
> Regards,
> Ariel
>
> 2018-09-10 17:29 GMT-03:00 MacQueen, Don <[hidden email]>:
>
> > If all of your data frames had enough points then this should work:
> >
> > myfun1 <- function(le) {
> >   list(one=geosphere::centroid(coordinates(le$one)),
> >        two=geosphere::centroid(coordinates(le$two))
> >        )
> > }
> >
> > lapply(ct, myfun1)
> >
> > To handle the one point case, try the following, which I think works with
> > your example data, but is untested if there just two points.
> > It also assumes the data frames are named 'one' and 'two', and will
> ignore
> > any others. To handle other names, I think you
> > could replace  $one  with  [[1]]  and  $two  with  [[2]]  .
> >
> > myfun <- function(le) {
> >   list(one=if (length(le$one)>1) geosphere::centroid(coordinates(le$one))
> > else coodinates(le$one),
> >        two=if (length(le$two)>1) geosphere::centroid(coordinates(le$two))
> > else coordinates(le$two)
> >        )
> > }
> >
> > lapply(ct, myfun)
> >
> >
> > To see a little bit of what's going on, try
> >
> >   myfun(ct[[1]])
> >   myfun1(ct[[1]])
> >   myfun1(ct[[2]])
> >   myfun(ct[[2]])
> >
> >
> >
> > I would also verify that the length method for objects of class
> > SpatialPoints does return the number of points, as it appears to:
> >
> > > class(ct$a$two)
> > [1] "SpatialPoints"
> > attr(,"package")
> > [1] "sp"
> >
> > > length(ct$a$one)
> > [1] 3
> > > length(ct$a$two)
> > [1] 3
> > > length(ct$a$two)
> > [1] 3
> > > length(ct$b$two)
> > [1] 1
> >
> > --
> > Don MacQueen
> > Lawrence Livermore National Laboratory
> > 7000 East Ave., L-627
> > Livermore, CA 94550
> > 925-423-1062
> > Lab cell 925-724-7509
> >
> >
> >
> > On 9/10/18, 10:56 AM, "R-sig-Geo on behalf of Ariel Fuentesdi" <
> > [hidden email] on behalf of [hidden email]>
> > wrote:
> >
> >     Hi everyone,
> >
> >     I have a list of coordinates called "ct" and I want to extract the
> >     centroids of each sublist, but it only works when it has only 3 or
> more
> >     points. In ct$b$two only has one point; in that case, I would rescue
> > that
> >     point as If it were a centroid.
> >
> >     This is what I did:
> >
> >     library(dplyr)
> >     library(geosphere)
> >
> >     ct <- list(a = list(one = data.frame(lon = c(-180, -160, -60), lat =
> > c(-20,
> >     5, 0)),
> >                         two = data.frame(lon = c(-18, -16, -6), lat =
> > c(-2, 50,
> >     10))),
> >                b = list(one = data.frame(lon = c(-9, -8, -3), lat = c(-1,
> > 25,
> >     5)),
> >                         two = data.frame(lon = c(-90), lat = c(-1))))
> >
> >     coordinates(ct$a$one) <- ~lon+lat
> >     coordinates(ct$a$two) <- ~lon+lat
> >     coordinates(ct$b$one) <- ~lon+lat
> >     coordinates(ct$b$two) <- ~lon+lat
> >
> >     s <- 1:length(ct)
> >     ctply <- list()
> >     for (i in s){
> >      for (j in 1:length(ct[[i]])) {
> >       ctply[[i]][j] <- ifelse(test = lengths(ct[[i]][1]) > 2, yes =
> > sapply(X =
> >     ct[[i]][j],
> >         FUN = function(y) geosphere::centroid(slot(y, "coords"))), no =
> >     ct[[i]][j] )
> >      }
> >     }
> >
> >     Thanks in advance
> >
> >     Regards,
> >     Ariel
> >
> >         [[alternative HTML version deleted]]
> >
> >     _______________________________________________
> >     R-sig-Geo mailing list
> >     [hidden email]
> >     https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
> >
> >
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

--
Vijay Lulla

Assistant Professor,
Dept. of Geography, IUPUI
425 University Blvd, CA-207C.
Indianapolis, IN-46202
[hidden email]
ORCID: https://orcid.org/0000-0002-0823-2522
Webpage: http://vijaylulla.com

        [[alternative HTML version deleted]]

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

Re: help: Problem getting centroids inside lists

Wed, 09/12/2018 - 13:41
I went to a more general approach, that is:

ctply <- lapply(X = s, FUN = function(x) sapply(X = ct[[x]],
          FUN = function(y) if(length(y)>1) geosphere::centroid(slot(y,
"coords"))
          else sp::coordinates(slot(y, "coords"))))

But now I want to add the case when they are only two elements. The dataset
will be:

ct <- list(a = list(one = data.frame(lon = c(-180, -160, -60), lat = c(-20,
5, 0)),
                    two = data.frame(lon = c(-18, -16), lat = c(-2, 50))),
           b = list(one = data.frame(lon = c(-9, -8, -3), lat = c(-1, 25,
5)),
                    two = data.frame(lon = c(-90), lat = c(-1))))
coordinates(ct$a$one) <- ~lon+lat
coordinates(ct$a$two) <- ~lon+lat
coordinates(ct$b$one) <- ~lon+lat
coordinates(ct$b$two) <- ~lon+lat

And I did the following but it doesn't work:

ctply <- lapply(X = s, FUN = function(x) sapply(X = ct[[x]],
                                                FUN = function(y)
if(length(y)>2) geosphere::centroid(slot(y, "coords"))
                                                else if (length(y) == 1)
sp::coordinates(slot(y, "coords"))
                                                else
mean(sp::coordinates(slot(y, "coords")))))

I need the result of each element of the list will be a matrix of two rows
per column (named: one, two). How do I fix it?

Regards,
Ariel

2018-09-10 17:29 GMT-03:00 MacQueen, Don <[hidden email]>:

> If all of your data frames had enough points then this should work:
>
> myfun1 <- function(le) {
>   list(one=geosphere::centroid(coordinates(le$one)),
>        two=geosphere::centroid(coordinates(le$two))
>        )
> }
>
> lapply(ct, myfun1)
>
> To handle the one point case, try the following, which I think works with
> your example data, but is untested if there just two points.
> It also assumes the data frames are named 'one' and 'two', and will ignore
> any others. To handle other names, I think you
> could replace  $one  with  [[1]]  and  $two  with  [[2]]  .
>
> myfun <- function(le) {
>   list(one=if (length(le$one)>1) geosphere::centroid(coordinates(le$one))
> else coodinates(le$one),
>        two=if (length(le$two)>1) geosphere::centroid(coordinates(le$two))
> else coordinates(le$two)
>        )
> }
>
> lapply(ct, myfun)
>
>
> To see a little bit of what's going on, try
>
>   myfun(ct[[1]])
>   myfun1(ct[[1]])
>   myfun1(ct[[2]])
>   myfun(ct[[2]])
>
>
>
> I would also verify that the length method for objects of class
> SpatialPoints does return the number of points, as it appears to:
>
> > class(ct$a$two)
> [1] "SpatialPoints"
> attr(,"package")
> [1] "sp"
>
> > length(ct$a$one)
> [1] 3
> > length(ct$a$two)
> [1] 3
> > length(ct$a$two)
> [1] 3
> > length(ct$b$two)
> [1] 1
>
> --
> Don MacQueen
> Lawrence Livermore National Laboratory
> 7000 East Ave., L-627
> Livermore, CA 94550
> 925-423-1062
> Lab cell 925-724-7509
>
>
>
> On 9/10/18, 10:56 AM, "R-sig-Geo on behalf of Ariel Fuentesdi" <
> [hidden email] on behalf of [hidden email]>
> wrote:
>
>     Hi everyone,
>
>     I have a list of coordinates called "ct" and I want to extract the
>     centroids of each sublist, but it only works when it has only 3 or more
>     points. In ct$b$two only has one point; in that case, I would rescue
> that
>     point as If it were a centroid.
>
>     This is what I did:
>
>     library(dplyr)
>     library(geosphere)
>
>     ct <- list(a = list(one = data.frame(lon = c(-180, -160, -60), lat =
> c(-20,
>     5, 0)),
>                         two = data.frame(lon = c(-18, -16, -6), lat =
> c(-2, 50,
>     10))),
>                b = list(one = data.frame(lon = c(-9, -8, -3), lat = c(-1,
> 25,
>     5)),
>                         two = data.frame(lon = c(-90), lat = c(-1))))
>
>     coordinates(ct$a$one) <- ~lon+lat
>     coordinates(ct$a$two) <- ~lon+lat
>     coordinates(ct$b$one) <- ~lon+lat
>     coordinates(ct$b$two) <- ~lon+lat
>
>     s <- 1:length(ct)
>     ctply <- list()
>     for (i in s){
>      for (j in 1:length(ct[[i]])) {
>       ctply[[i]][j] <- ifelse(test = lengths(ct[[i]][1]) > 2, yes =
> sapply(X =
>     ct[[i]][j],
>         FUN = function(y) geosphere::centroid(slot(y, "coords"))), no =
>     ct[[i]][j] )
>      }
>     }
>
>     Thanks in advance
>
>     Regards,
>     Ariel
>
>         [[alternative HTML version deleted]]
>
>     _______________________________________________
>     R-sig-Geo mailing list
>     [hidden email]
>     https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
>
        [[alternative HTML version deleted]]

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

Help with simple Map of US states to predefined regions

Wed, 09/12/2018 - 13:27
Hi

I have this df with three columns ProviderState, ProviderStateCode, ProviderRegion I wanted to use to create a simple 5 color map

I have reviewed fiftystater pkg and map pkg but not sure how to simply take these three columns and plot a simple 5 color map based on the Region the state is in?

After looking at these and trying to apply these ideas to my data
https://cran.r-project.org/web/packages/fiftystater/vignettes/fiftystater.html
https://cran.r-project.org/web/packages/maps/maps.pdf

I found tutorial at: https://uchicagoconsulting.wordpress.com/tag/r-ggplot2-maps-visualization/

I used the tutorial data and subset in my regions

So now I have come up with the 5 segmented maps and my question becomes how to put this all into one map of the US?


install.packages("maps")
library(maps)
library(ggplot2)

#load us map data
all_states <- map_data("state") View(all_states)
#plot all states with ggplot
p <- ggplot()
p <- p + geom_polygon( data=all_states, aes(x=long, y=lat, group = group),colour="white", fill="blue" )
p

#http://sape.inf.usi.ch/quick-reference/ggplot2/colour

#Pacificstates
Pacificstates <- subset(all_states, region %in% c( "alaska", "arizona", "california", "hawaii", "nevada", "oregon","washington") )
p <- ggplot()
p <- p + geom_polygon( data=Pacificstates, aes(x=long, y=lat, group = group),colour="white", fill="deepskyblue4" ) +
  labs(title = "Pacificstates")
p

#Frontierstates
Frontierstates <- subset(all_states, region %in% c( "colorado", "idaho", "kansas", "montana", "new mexico", "oklahoma","texas", "utah", "wyoming") )
p <- ggplot()
p <- p + geom_polygon( data=Frontierstates, aes(x=long, y=lat, group = group),colour="white", fill="dodgerblue1" ) +
      labs(title = "FrontierStates")
p

#Midweststates
Midweststates <- subset(all_states, region %in% c( "iowa", "illinois", "indiana", "michigan", "minnesota", "missouri","north dakota", "nebraska", "ohio","south dakota","wisconsin") )
p <- ggplot()
p <- p + geom_polygon( data=Midweststates, aes(x=long, y=lat, group = group),colour="white", fill="dodgerblue1" ) +
      labs(title = "MidwestStates")
p

#Southernstates
Southernstates <- subset(all_states, region %in% c( "alabama", "arkansas", "florida", "georgia", "kentucky", "louisiana","mississippi"
                                                    ,"north carolina", "south carolina","tennessee","virginia","west virginia") )
p <- ggplot()
p <- p + geom_polygon( data=Southernstates, aes(x=long, y=lat, group = group),colour="white", fill="royalblue2" ) +
  labs(title = "Southernstates")
p

# Northeaststates
Northeaststates <- subset(all_states, region %in% c( "connecticut", "district of columbia", "delaware", "massachusetts", "maryland", "maine","new hampshire"
                                                     , "new jersey", "new york","pennsylvania","rhode island","vermont") )
p <- ggplot()
p <- p + geom_polygon( data=Northeaststates, aes(x=long, y=lat, group = group),colour="white", fill="dodgerblue4" ) +
  labs(title = "Northeaststates")
p


#here is my my data but not used above

str(Map1)
Classes 'tbl_df', 'tbl' and 'data.frame':    54 obs. of  3 variables:
$ ProviderState    : chr  "ALASKA" "ALABAMA" "ARKANSAS" "ARIZONA" ...
$ ProviderStateCode: chr  "AK" "AL" "AR" "AZ" ...
$ ProviderRegion   : chr  "Pacific" "South" "South" "Pacific" ...
- attr(*, "spec")=List of 2
  ..$ cols   :List of 3
  .. ..$ ProviderState    : list()
  .. .. ..- attr(*, "class")= chr  "collector_character" "collector"
  .. ..$ ProviderStateCode: list()
  .. .. ..- attr(*, "class")= chr  "collector_character" "collector"
  .. ..$ ProviderRegion   : list()
  .. .. ..- attr(*, "class")= chr  "collector_character" "collector"
  ..$ default: list()
  .. ..- attr(*, "class")= chr  "collector_guess" "collector"
  ..- attr(*, "class")= chr "col_spec"

dput(Map1)
structure(list(ProviderState = c("ALASKA", "ALABAMA", "ARKANSAS",
"ARIZONA", "CALIFORNIA", "COLORADO", "CONNECTICUT", "DISTRICT OF COLUMBIA",
"DELAWARE", "FLORIDA", "GEORGIA", "GUAM", "HAWAII", "IOWA", "IDAHO",
"ILLINOIS", "INDIANA", "KANSAS", "KENTUCKY", "LOUISIANA", "MASSACHUSETTS",
"MARYLAND", "MAINE", "MICHIGAN", "MINNESOTA", "MISSOURI", "MISSISSIPPI",
"MONTANA", "NORTH CAROLINA", "NORTH DAKOTA", "NEBRASKA", "NEW HAMPSHIRE",
"NEW JERSEY", "NEW MEXICO", "NEVADA", "NEW YORK", "OHIO", "OKLAHOMA",
"OREGON", "PENNSYLVANIA", "PUERTO RICO", "RHODE ISLAND", "SOUTH CAROLINA",
"SOUTH DAKOTA", "TENNESSEE", "TEXAS", "UTAH", "VIRGINIA", "VIRGIN ISLANDS",
"VERMONT", "WASHINGTON", "WISCONSIN", "WEST VIRGINIA", "WYOMING"
), ProviderStateCode = c("AK", "AL", "AR", "AZ", "CA", "CO",
"CT", "DC", "DE", "FL", "GA", "GU", "HI", "IA", "ID", "IL", "IN",
"KS", "KY", "LA", "MA", "MD", "ME", "MI", "MN", "MO", "MS", "MT",
"NC", "ND", "NE", "NH", "NJ", "NM", "NV", "NY", "OH", "OK", "OR",
"PA", "PR", "RI", "SC", "SD", "TN", "TX", "UT", "VA", "VI", "VT",
"WA", "WI", "WV", "WY"), ProviderRegion = c("Pacific", "South",
"South", "Pacific", "Pacific", "Frontier", "Northeast", "Northeast",
"Northeast", "South", "South", "Pacific", "Pacific", "Midwest",
"Frontier", "Midwest", "Midwest", "Frontier", "South", "South",
"Northeast", "Northeast", "Northeast", "Midwest", "Midwest",
"Midwest", "South", "Frontier", "South", "Midwest", "Midwest",
"Northeast", "Northeast", "Frontier", "Pacific", "Northeast",
"Midwest", "Frontier", "Pacific", "Northeast", "Northeast", "Northeast",
"South", "Midwest", "South", "Frontier", "Frontier", "South",
"Northeast", "Northeast", "Pacific", "Midwest", "South", "Frontier"
)), row.names = c(NA, -54L), class = c("tbl_df", "tbl", "data.frame"
), spec = structure(list(cols = list(ProviderState = structure(list(), class = c("collector_character",
"collector")), ProviderStateCode = structure(list(), class = c("collector_character",
"collector")), ProviderRegion = structure(list(), class = c("collector_character",
"collector"))), default = structure(list(), class = c("collector_guess",
"collector"))), class = "col_spec"))

Thank you for any suggestions.

WHP






Confidentiality Notice This message is sent from Zelis. ...{{dropped:15}}

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

Re: Leaflet map nested in RShiny App - Improving speed & portability

Wed, 09/12/2018 - 13:26
Hey Erin,

I agree with everybody who answered, but just for completeness and FYI:

Another approach for sharing a Shiny app offline as a (Windows) Desktop
app might be using Shiny with "R portable".

For further details please check for example:
https://www.r-bloggers.com/deploying-desktop-apps-with-r/

Cheers!

Harry

Ps. If anybody thinks this is a very bad idea, please let me know,
because I'm planing to develop something like that ... Thanks!


Am 05.09.2018 um 01:56 schrieb Erin Stearns:
> Hello all!
>
> I hope this message finds you all well!
>
> I have 2 questions pertaining to the creation of interactive maps via
> Leaflet nested inside an RShiny app. One question has to do with
> computation while the other has to do with sharing/off-line interactivity.
>
> *Project description:*
> I am creating a global map depicting binary malaria risk (at risk, not at
> risk) at the Admin 2 level(current state only uses 5 countries and can be
> found here <https://erstearns.shinyapps.io/malariarisk5/>).  I am using an
> ESRI base map, then a polygons shapefile containing geometry and attributes
> (geographical hierarchy & risk).
>
> *Computation question*
> As you see, the RShiny app takes quite a bit of time to render. Does anyone
> have any suggestions for improving this? As previously said, this version
> only contains 5 countries, thus I cannot continue with my current method to
> reach a global map. I have considered finding centroids of all Admin 2
> polygons and retaining attribute information here, then rasterizing the
> malaria risk shapefile for visualization and using the 2 instead of a
> single shapefile with polygon boundaries and attributes.
>
> *Sharing the app/offline interactivity*
> I am planning to share this with people who likely do not have R installed
> on their laptops nor have they ever coded. Does anyone have any suggestions
> for the best way to do this while retaining interactivity?
>
> Thank you all, any insight is greatly appreciated.
>
> Best,
> Erin
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: subsetting a spatial polygons

Wed, 09/12/2018 - 08:08
Thanks Michael and Tim, thanks for the attention

The first thing I did was I try the subsets using subset.

subset(polys,
coordinates(polys)[,1]>=-46 & coordinates(polys)[,1]<=-45 &
coordinates(polys)[,2]>=-25 & coordinates(polys)[,2]<=-24)

But it seems to work only with Spatial*DataFrame objects. Then I (unluckily)
gave up the x-y approach.

With Mike's advice I tried again as

plot(
    polys[
        coordinates(polys)[,1]>=-46 & coordinates(polys)[,1]<=-45 &
        coordinates(polys)[,2]>=-25 & coordinates(polys)[,2]<=-24,]
,add=T,col="green")

and I got what I wanted in a more "elegant" fashion.

Once more thanks a lot for the attention.

--
Antônio Olinto Ávila da Silva
Instituto de Pesca (Fisheries Institute)
São Paulo, Brasil

Em qua, 12 de set de 2018 às 05:05, Michael Sumner <[hidden email]>
escreveu:

> crop doesn't select it actually cuts on the boundary, for simple shapes I
> use coordinates (for centroids in easy matrix) and select with [ using
> tests on x and y.
>
> Cheers, Mike
>
> On Wed, Sep 12, 2018, 15:29 Tim Appelhans <[hidden email]> wrote:
>
> > Antonio,
> >
> > I am not sure why you think that your solution is not very elegant.
> > In case you want to have more visual control over the subsetting, you
> > could try mapedit:
> >
> > library(mapedit)
> > myselection = selectFeatures(polys, mode = "draw")
> >
> > which will let you draw a e.g. rectangle and only return those features
> > that intersect it.
> >
> > Best
> > Tim
> >
> >
> > On 09/12/2018 04:18 AM, Antonio Silva wrote:
> > > Dear list users
> > >
> > > I have a SpatialPolygons with several squares. How to subset it to have
> > > only the squares between given latitudes and longitudes?
> > >
> > > In the example
> > >
> > > library(sp)
> > > library(raster)
> > > grd <-
> > >
> >
> GridTopology(cellcentre.offset=c(-47.75,-25.416667),cellsize=c(10/60,10/60),cells.dim=c(23,12))
> > > polys <- as.SpatialPolygons.GridTopology(grd)
> > > proj4string(polys) <- CRS("+proj=longlat +datum=WGS84")
> > > plot(polys,axes=T)
> > >
> > > How to select only the squares, let's say, between 24-25°S and 45-46°W?
> > >
> > > The farthest I went was:
> > >
> > > e <- extent(-45.9,-45.1,-24.9,-24.1) # which is not very elegant
> > > mask <- crop(polys,e)
> > > polys2 <- polys[mask,]
> > > plot(polys2,add=T,col="green")
> > >
> > > Thanks a lot. Best regards
> > >
> >
> > _______________________________________________
> > 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
>
        [[alternative HTML version deleted]]

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

Re: subsetting a spatial polygons

Wed, 09/12/2018 - 03:04
crop doesn't select it actually cuts on the boundary, for simple shapes I
use coordinates (for centroids in easy matrix) and select with [ using
tests on x and y.

Cheers, Mike

On Wed, Sep 12, 2018, 15:29 Tim Appelhans <[hidden email]> wrote:

> Antonio,
>
> I am not sure why you think that your solution is not very elegant.
> In case you want to have more visual control over the subsetting, you
> could try mapedit:
>
> library(mapedit)
> myselection = selectFeatures(polys, mode = "draw")
>
> which will let you draw a e.g. rectangle and only return those features
> that intersect it.
>
> Best
> Tim
>
>
> On 09/12/2018 04:18 AM, Antonio Silva wrote:
> > Dear list users
> >
> > I have a SpatialPolygons with several squares. How to subset it to have
> > only the squares between given latitudes and longitudes?
> >
> > In the example
> >
> > library(sp)
> > library(raster)
> > grd <-
> >
> GridTopology(cellcentre.offset=c(-47.75,-25.416667),cellsize=c(10/60,10/60),cells.dim=c(23,12))
> > polys <- as.SpatialPolygons.GridTopology(grd)
> > proj4string(polys) <- CRS("+proj=longlat +datum=WGS84")
> > plot(polys,axes=T)
> >
> > How to select only the squares, let's say, between 24-25°S and 45-46°W?
> >
> > The farthest I went was:
> >
> > e <- extent(-45.9,-45.1,-24.9,-24.1) # which is not very elegant
> > mask <- crop(polys,e)
> > polys2 <- polys[mask,]
> > plot(polys2,add=T,col="green")
> >
> > Thanks a lot. Best regards
> >
>
> _______________________________________________
> 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

Re: subsetting a spatial polygons

Wed, 09/12/2018 - 00:28
Antonio,

I am not sure why you think that your solution is not very elegant.
In case you want to have more visual control over the subsetting, you
could try mapedit:

library(mapedit)
myselection = selectFeatures(polys, mode = "draw")

which will let you draw a e.g. rectangle and only return those features
that intersect it.

Best
Tim


On 09/12/2018 04:18 AM, Antonio Silva wrote:
> Dear list users
>
> I have a SpatialPolygons with several squares. How to subset it to have
> only the squares between given latitudes and longitudes?
>
> In the example
>
> library(sp)
> library(raster)
> grd <-
> GridTopology(cellcentre.offset=c(-47.75,-25.416667),cellsize=c(10/60,10/60),cells.dim=c(23,12))
> polys <- as.SpatialPolygons.GridTopology(grd)
> proj4string(polys) <- CRS("+proj=longlat +datum=WGS84")
> plot(polys,axes=T)
>
> How to select only the squares, let's say, between 24-25°S and 45-46°W?
>
> The farthest I went was:
>
> e <- extent(-45.9,-45.1,-24.9,-24.1) # which is not very elegant
> mask <- crop(polys,e)
> polys2 <- polys[mask,]
> plot(polys2,add=T,col="green")
>
> Thanks a lot. Best regards
>
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

subsetting a spatial polygons

Tue, 09/11/2018 - 21:18
Dear list users

I have a SpatialPolygons with several squares. How to subset it to have
only the squares between given latitudes and longitudes?

In the example

library(sp)
library(raster)
grd <-
GridTopology(cellcentre.offset=c(-47.75,-25.416667),cellsize=c(10/60,10/60),cells.dim=c(23,12))
polys <- as.SpatialPolygons.GridTopology(grd)
proj4string(polys) <- CRS("+proj=longlat +datum=WGS84")
plot(polys,axes=T)

How to select only the squares, let's say, between 24-25°S and 45-46°W?

The farthest I went was:

e <- extent(-45.9,-45.1,-24.9,-24.1) # which is not very elegant
mask <- crop(polys,e)
polys2 <- polys[mask,]
plot(polys2,add=T,col="green")

Thanks a lot. Best regards

--
Antônio Olinto Ávila da Silva
Instituto de Pesca (Fisheries Institute)
São Paulo, Brasil

        [[alternative HTML version deleted]]

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

Re: help: Problem getting centroids inside lists

Mon, 09/10/2018 - 15:29
If all of your data frames had enough points then this should work:

myfun1 <- function(le) {
  list(one=geosphere::centroid(coordinates(le$one)),
       two=geosphere::centroid(coordinates(le$two))
       )
}

lapply(ct, myfun1)

To handle the one point case, try the following, which I think works with your example data, but is untested if there just two points.
It also assumes the data frames are named 'one' and 'two', and will ignore any others. To handle other names, I think you
could replace  $one  with  [[1]]  and  $two  with  [[2]]  .

myfun <- function(le) {
  list(one=if (length(le$one)>1) geosphere::centroid(coordinates(le$one)) else coodinates(le$one),
       two=if (length(le$two)>1) geosphere::centroid(coordinates(le$two)) else coordinates(le$two)
       )
}

lapply(ct, myfun)


To see a little bit of what's going on, try

  myfun(ct[[1]])
  myfun1(ct[[1]])
  myfun1(ct[[2]])
  myfun(ct[[2]])



I would also verify that the length method for objects of class SpatialPoints does return the number of points, as it appears to:

> class(ct$a$two)
[1] "SpatialPoints"
attr(,"package")
[1] "sp"

> length(ct$a$one)
[1] 3
> length(ct$a$two)
[1] 3
> length(ct$a$two)
[1] 3
> length(ct$b$two)
[1] 1

--
Don MacQueen
Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062
Lab cell 925-724-7509
 
 

On 9/10/18, 10:56 AM, "R-sig-Geo on behalf of Ariel Fuentesdi" <[hidden email] on behalf of [hidden email]> wrote:

    Hi everyone,
   
    I have a list of coordinates called "ct" and I want to extract the
    centroids of each sublist, but it only works when it has only 3 or more
    points. In ct$b$two only has one point; in that case, I would rescue that
    point as If it were a centroid.
   
    This is what I did:
   
    library(dplyr)
    library(geosphere)
   
    ct <- list(a = list(one = data.frame(lon = c(-180, -160, -60), lat = c(-20,
    5, 0)),
                        two = data.frame(lon = c(-18, -16, -6), lat = c(-2, 50,
    10))),
               b = list(one = data.frame(lon = c(-9, -8, -3), lat = c(-1, 25,
    5)),
                        two = data.frame(lon = c(-90), lat = c(-1))))
   
    coordinates(ct$a$one) <- ~lon+lat
    coordinates(ct$a$two) <- ~lon+lat
    coordinates(ct$b$one) <- ~lon+lat
    coordinates(ct$b$two) <- ~lon+lat
   
    s <- 1:length(ct)
    ctply <- list()
    for (i in s){
     for (j in 1:length(ct[[i]])) {
      ctply[[i]][j] <- ifelse(test = lengths(ct[[i]][1]) > 2, yes = sapply(X =
    ct[[i]][j],
        FUN = function(y) geosphere::centroid(slot(y, "coords"))), no =
    ct[[i]][j] )
     }
    }
   
    Thanks in advance
   
    Regards,
    Ariel
   
    [[alternative HTML version deleted]]
   
    _______________________________________________
    R-sig-Geo mailing list
    [hidden email]
    https://stat.ethz.ch/mailman/listinfo/r-sig-geo
   

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

help: Problem getting centroids inside lists

Mon, 09/10/2018 - 12:56
Hi everyone,

I have a list of coordinates called "ct" and I want to extract the
centroids of each sublist, but it only works when it has only 3 or more
points. In ct$b$two only has one point; in that case, I would rescue that
point as If it were a centroid.

This is what I did:

library(dplyr)
library(geosphere)

ct <- list(a = list(one = data.frame(lon = c(-180, -160, -60), lat = c(-20,
5, 0)),
                    two = data.frame(lon = c(-18, -16, -6), lat = c(-2, 50,
10))),
           b = list(one = data.frame(lon = c(-9, -8, -3), lat = c(-1, 25,
5)),
                    two = data.frame(lon = c(-90), lat = c(-1))))

coordinates(ct$a$one) <- ~lon+lat
coordinates(ct$a$two) <- ~lon+lat
coordinates(ct$b$one) <- ~lon+lat
coordinates(ct$b$two) <- ~lon+lat

s <- 1:length(ct)
ctply <- list()
for (i in s){
 for (j in 1:length(ct[[i]])) {
  ctply[[i]][j] <- ifelse(test = lengths(ct[[i]][1]) > 2, yes = sapply(X =
ct[[i]][j],
    FUN = function(y) geosphere::centroid(slot(y, "coords"))), no =
ct[[i]][j] )
 }
}

Thanks in advance

Regards,
Ariel

        [[alternative HTML version deleted]]

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

Isohyet maps where the inverse distance method does not work correctly

Mon, 09/10/2018 - 07:48
Dear list users,
I implemented a code for creation of snowfall isohyets through the inverse distance method applied to a Digital Elevation Model.
Within a big area (approximatively 2.000 km2) I use only six stations, I would expect the snowfall esteem getting smaller with higher distances from the stations.
In my case the snowfall values remain constant all over the area.

This is the relevant part of my code (df_prec is the data frame that contains the coordinates of the stations and the snowfall cumulates).
Could somebody be so patient to spot where is my (big) mistake? (Here attached - as exception - the file picture1.pdf is an example of the result.)

Thank you for your attention and your help
Stefano Sofia



## LOADING SHAPEFILE FOR AREA5
Shape_area <- readOGR(dsn="area.shp", layer="area")
Shape_area_projection <- "+init=epsg:32633 +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
Shape_area <- spTransform(Shape_area, CRS(Shape_area_projection))

## LOADING DEM
area_DEM raster("dem.tif")
area_DEM <- projectRaster(area_DEM, crs = CRS("+init=epsg:32633"))

## EXTRACTING THE ELEVATION VALUES TO MY POINTS
df_prec$ExtractedElevationValues <- extract(x=area_DEM, y=df_prec)

## CREATING A NEW GRID
newgrid <- as(area_DEM, "SpatialGridDataFrame")
newgrid <- raster(newgrid)
names(newgrid) <- "ExtractedElevationValues"

## Inverse Distance Weighted
if (sum(df_prec$Cumulata) > 0)
{
  OK_rain <- gstat(formula=snowfall_cumulate~ExtractedElevationValues, data=df_prec, locations=newgrid)
  rain_rast <- interpolate(newgrid, OK_rain, xyOnly=FALSE, na.rm=FALSE)

} else
{
  rain_rast <- raster(df_prec)*0
}




         (oo)
--oOO--( )--OOo----------------
Stefano Sofia PhD
Area Meteorologica e  Area nivologica - Centro Funzionale
Servizio Protezione Civile - Regione Marche
Via del Colle Ameno 5
60126 Torrette di Ancona, Ancona
Uff: 071 806 7743
E-mail: [hidden email]
---Oo---------oO----------------

AVVISO IMPORTANTE: Questo messaggio di posta elettronica può contenere informazioni confidenziali, pertanto è destinato solo a persone autorizzate alla ricezione. I messaggi di posta elettronica per i client di Regione Marche possono contenere informazioni confidenziali e con privilegi legali. Se non si è il destinatario specificato, non leggere, copiare, inoltrare o archiviare questo messaggio. Se si è ricevuto questo messaggio per errore, inoltrarlo al mittente ed eliminarlo completamente dal sistema del proprio computer. Ai sensi dell’art. 6 della DGR n. 1394/2008 si segnala che, in caso di necessità ed urgenza, la risposta al presente messaggio di posta elettronica può essere visionata da persone estranee al destinatario.
IMPORTANT NOTICE: This e-mail message is intended to be received only by persons entitled to receive the confidential information it may contain. E-mail messages to clients of Regione Marche may contain information that is confidential and legally privileged. Please do not read, copy, forward, or store this message unless you are an intended recipient of it. If you have received this message in error, please forward it to the sender and delete it completely from your computer system.

--
Questo messaggio è stato analizzato con Libra ESVA ed è risultato non infetto.
This message has been checked by Libra ESVA and is believed to be clean.


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

picture1.pdf (22K) Download Attachment

Re: Leaflet map nested in RShiny App - Improving speed & portability

Thu, 09/06/2018 - 12:40
Thank you, Mike! Yes, I have used htmlwidgets, thank you for that! The
issue is the file begins unwieldy given the amount of data contained in my
leaflet app.

Yeah, it seems like there are a number of options, it's just trying to
determine the best is the tricky part.

Thanks again for the ideas - very much appreciated.

Best,
Erin



On Wed, Sep 5, 2018 at 11:02 AM Michael Treglia <[hidden email]> wrote:

> I'll just second Barry's idea in particular, to set up as a standalone
> webpage. You could even use QGIS and the QGIS2Web Plugin to create that,
> and host via GitHub pages or similar.
>
> From R, after creating a map via leaflet and similar packages, you can use
> htmlwidgets::saveWidget() to export as a standalone .html file if I recall
> correctly.
>
> The one thing regarding a standalone webpage is that if you have a lot of
> objects (especially complex ones), that can be a lot for a browser to
> handle (given the data are part of the html file).  Might be worth some
> quick experimentation, and simplifying polygons would help. (You could
> always create a quick landing page, even generated via rMarkdown, and
> having a link for maps by different regions or countries - then you could
> have a folder of .html files you could distribute, and users could just
> open the landing page, and navigate from there).
>
> Just some quick thoughts... Hope this helps.
> Mike T
>
> On Wed, Sep 5, 2018 at 1:17 PM Erin Stearns <[hidden email]> wrote:
>
>> Hello all,
>>
>> Thank you all very much for the great insight!
>>
>> *McCrea *- thank you very much, I will test using a geojson first, then
>> test after reducing geometry.
>>
>> *Tim* - thank you for the great breakdown and recommended priority list.
>> Ideally, I would like to be able to share the interactive map with
>> teammates as a file or something akin to it such that they can simply open
>> it and interact with the map. RInno is a great option, however I run a
>> linux machine, so will look into further, but may need to find another
>> option.
>>
>> *Roman* - the app is currently deployed to shinyapps.io. Thank you for
>> sharing about ShinyProxy -- so would this method require 1. Internet and
>> 2.
>> local installation (vs internal server)?
>>
>> *Barry* - wow, thank you for your response! Sounds like this would be the
>> best way to solve both issues. I am not as fluent with HTML and JS, but as
>> you say, there are likely great guides available to take this route.
>>
>> Thank you all again, this has been hugely helpful. I wish you all the best
>> and hope I can be of help to you at some point!
>>
>> Best,
>> Erin
>>
>>
>> On Wed, Sep 5, 2018 at 12:48 AM Barry Rowlingson <
>> [hidden email]> wrote:
>>
>> >
>> >
>> > On Wed, Sep 5, 2018 at 12:56 AM, Erin Stearns <[hidden email]>
>> > wrote:
>> >
>> >> Hello all!
>> >>
>> >> I hope this message finds you all well!
>> >>
>> >> I have 2 questions pertaining to the creation of interactive maps via
>> >> Leaflet nested inside an RShiny app. One question has to do with
>> >> computation while the other has to do with sharing/off-line
>> interactivity.
>> >>
>> >> *Computation question*
>> >> As you see, the RShiny app takes quite a bit of time to render. Does
>> >> anyone
>> >> have any suggestions for improving this? As previously said, this
>> version
>> >> only contains 5 countries, thus I cannot continue with my current
>> method
>> >> to
>> >> reach a global map. I have considered finding centroids of all Admin 2
>> >> polygons and retaining attribute information here, then rasterizing the
>> >> malaria risk shapefile for visualization and using the 2 instead of a
>> >> single shapefile with polygon boundaries and attributes.
>> >>
>> >>
>> > Unless you plan to add any computational functions to this map then I'd
>> > strongly recommend creating it as a standalone web app and not a shiny
>> app.
>> > This will enable you to use lots of useful Leaflet plugins for speeding
>> > things up, such as only showing country outlines at low zoom levels, and
>> > showing subdivisions only at high zoom levels. This *might* be possible
>> > with R's various leaflet packages but I'd go for full javascript
>> control.
>> >
>> > A standalone map would take its data from a JSON file or similar, and
>> you
>> > would then be writing R code that generated that. The mapping app
>> itself is
>> > written in HTML and JS with CSS styling. There are plenty of guides to
>> > web-based interactive mapping, starting with Leaflet.
>> >
>> >
>> >> *Sharing the app/offline interactivity*
>> >> I am planning to share this with people who likely do not have R
>> installed
>> >> on their laptops nor have they ever coded. Does anyone have any
>> >> suggestions
>> >> for the best way to do this while retaining interactivity?
>> >>
>> >>  Here's the big win of creating a standalone web map. You only have to
>> > distribute the HTML/CSS/JS and they can be viewed directly (or you also
>> > supply a tiny server that runs locally and only has to feed the files
>> on a
>> > localhost port). No need to have a shiny server anywhere, or install R.
>> Its
>> > simple and clean. It also needs no network connectivity, but you'll not
>> get
>> > a base map - but you could include a low or medium resolution basemap
>> > raster in your package.
>> >
>> > The only reason to need Shiny here would be if you wanted people to do
>> > something computational, like click on a bunch of polygons and then fit
>> a
>> > linear model to the selection, since that would require a round-trip to
>> the
>> > server for R to compute the fit. (although I suspect there's a JS
>> package
>> > for linear modelling.... you can do ML in JS these days...)
>> >
>> >
>> >
>> >> Thank you all, any insight is greatly appreciated.
>> >>
>> >> Best,
>> >> Erin
>> >>
>> >>         [[alternative HTML version deleted]]
>> >>
>> >> _______________________________________________
>> >> R-sig-Geo mailing list
>> >> [hidden email]
>> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >>
>> >
>> >
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
        [[alternative HTML version deleted]]

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

Re: Leaflet map nested in RShiny App - Improving speed & portability

Wed, 09/05/2018 - 13:02
I'll just second Barry's idea in particular, to set up as a standalone
webpage. You could even use QGIS and the QGIS2Web Plugin to create that,
and host via GitHub pages or similar.

From R, after creating a map via leaflet and similar packages, you can use
htmlwidgets::saveWidget() to export as a standalone .html file if I recall
correctly.

The one thing regarding a standalone webpage is that if you have a lot of
objects (especially complex ones), that can be a lot for a browser to
handle (given the data are part of the html file).  Might be worth some
quick experimentation, and simplifying polygons would help. (You could
always create a quick landing page, even generated via rMarkdown, and
having a link for maps by different regions or countries - then you could
have a folder of .html files you could distribute, and users could just
open the landing page, and navigate from there).

Just some quick thoughts... Hope this helps.
Mike T

On Wed, Sep 5, 2018 at 1:17 PM Erin Stearns <[hidden email]> wrote:

> Hello all,
>
> Thank you all very much for the great insight!
>
> *McCrea *- thank you very much, I will test using a geojson first, then
> test after reducing geometry.
>
> *Tim* - thank you for the great breakdown and recommended priority list.
> Ideally, I would like to be able to share the interactive map with
> teammates as a file or something akin to it such that they can simply open
> it and interact with the map. RInno is a great option, however I run a
> linux machine, so will look into further, but may need to find another
> option.
>
> *Roman* - the app is currently deployed to shinyapps.io. Thank you for
> sharing about ShinyProxy -- so would this method require 1. Internet and 2.
> local installation (vs internal server)?
>
> *Barry* - wow, thank you for your response! Sounds like this would be the
> best way to solve both issues. I am not as fluent with HTML and JS, but as
> you say, there are likely great guides available to take this route.
>
> Thank you all again, this has been hugely helpful. I wish you all the best
> and hope I can be of help to you at some point!
>
> Best,
> Erin
>
>
> On Wed, Sep 5, 2018 at 12:48 AM Barry Rowlingson <
> [hidden email]> wrote:
>
> >
> >
> > On Wed, Sep 5, 2018 at 12:56 AM, Erin Stearns <[hidden email]>
> > wrote:
> >
> >> Hello all!
> >>
> >> I hope this message finds you all well!
> >>
> >> I have 2 questions pertaining to the creation of interactive maps via
> >> Leaflet nested inside an RShiny app. One question has to do with
> >> computation while the other has to do with sharing/off-line
> interactivity.
> >>
> >> *Computation question*
> >> As you see, the RShiny app takes quite a bit of time to render. Does
> >> anyone
> >> have any suggestions for improving this? As previously said, this
> version
> >> only contains 5 countries, thus I cannot continue with my current method
> >> to
> >> reach a global map. I have considered finding centroids of all Admin 2
> >> polygons and retaining attribute information here, then rasterizing the
> >> malaria risk shapefile for visualization and using the 2 instead of a
> >> single shapefile with polygon boundaries and attributes.
> >>
> >>
> > Unless you plan to add any computational functions to this map then I'd
> > strongly recommend creating it as a standalone web app and not a shiny
> app.
> > This will enable you to use lots of useful Leaflet plugins for speeding
> > things up, such as only showing country outlines at low zoom levels, and
> > showing subdivisions only at high zoom levels. This *might* be possible
> > with R's various leaflet packages but I'd go for full javascript control.
> >
> > A standalone map would take its data from a JSON file or similar, and you
> > would then be writing R code that generated that. The mapping app itself
> is
> > written in HTML and JS with CSS styling. There are plenty of guides to
> > web-based interactive mapping, starting with Leaflet.
> >
> >
> >> *Sharing the app/offline interactivity*
> >> I am planning to share this with people who likely do not have R
> installed
> >> on their laptops nor have they ever coded. Does anyone have any
> >> suggestions
> >> for the best way to do this while retaining interactivity?
> >>
> >>  Here's the big win of creating a standalone web map. You only have to
> > distribute the HTML/CSS/JS and they can be viewed directly (or you also
> > supply a tiny server that runs locally and only has to feed the files on
> a
> > localhost port). No need to have a shiny server anywhere, or install R.
> Its
> > simple and clean. It also needs no network connectivity, but you'll not
> get
> > a base map - but you could include a low or medium resolution basemap
> > raster in your package.
> >
> > The only reason to need Shiny here would be if you wanted people to do
> > something computational, like click on a bunch of polygons and then fit a
> > linear model to the selection, since that would require a round-trip to
> the
> > server for R to compute the fit. (although I suspect there's a JS package
> > for linear modelling.... you can do ML in JS these days...)
> >
> >
> >
> >> Thank you all, any insight is greatly appreciated.
> >>
> >> Best,
> >> Erin
> >>
> >>         [[alternative HTML version deleted]]
> >>
> >> _______________________________________________
> >> R-sig-Geo mailing list
> >> [hidden email]
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>
> >
> >
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
        [[alternative HTML version deleted]]

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

Re: Leaflet map nested in RShiny App - Improving speed & portability

Wed, 09/05/2018 - 12:17
Hello all,

Thank you all very much for the great insight!

*McCrea *- thank you very much, I will test using a geojson first, then
test after reducing geometry.

*Tim* - thank you for the great breakdown and recommended priority list.
Ideally, I would like to be able to share the interactive map with
teammates as a file or something akin to it such that they can simply open
it and interact with the map. RInno is a great option, however I run a
linux machine, so will look into further, but may need to find another
option.

*Roman* - the app is currently deployed to shinyapps.io. Thank you for
sharing about ShinyProxy -- so would this method require 1. Internet and 2.
local installation (vs internal server)?

*Barry* - wow, thank you for your response! Sounds like this would be the
best way to solve both issues. I am not as fluent with HTML and JS, but as
you say, there are likely great guides available to take this route.

Thank you all again, this has been hugely helpful. I wish you all the best
and hope I can be of help to you at some point!

Best,
Erin


On Wed, Sep 5, 2018 at 12:48 AM Barry Rowlingson <
[hidden email]> wrote:

>
>
> On Wed, Sep 5, 2018 at 12:56 AM, Erin Stearns <[hidden email]>
> wrote:
>
>> Hello all!
>>
>> I hope this message finds you all well!
>>
>> I have 2 questions pertaining to the creation of interactive maps via
>> Leaflet nested inside an RShiny app. One question has to do with
>> computation while the other has to do with sharing/off-line interactivity.
>>
>> *Computation question*
>> As you see, the RShiny app takes quite a bit of time to render. Does
>> anyone
>> have any suggestions for improving this? As previously said, this version
>> only contains 5 countries, thus I cannot continue with my current method
>> to
>> reach a global map. I have considered finding centroids of all Admin 2
>> polygons and retaining attribute information here, then rasterizing the
>> malaria risk shapefile for visualization and using the 2 instead of a
>> single shapefile with polygon boundaries and attributes.
>>
>>
> Unless you plan to add any computational functions to this map then I'd
> strongly recommend creating it as a standalone web app and not a shiny app.
> This will enable you to use lots of useful Leaflet plugins for speeding
> things up, such as only showing country outlines at low zoom levels, and
> showing subdivisions only at high zoom levels. This *might* be possible
> with R's various leaflet packages but I'd go for full javascript control.
>
> A standalone map would take its data from a JSON file or similar, and you
> would then be writing R code that generated that. The mapping app itself is
> written in HTML and JS with CSS styling. There are plenty of guides to
> web-based interactive mapping, starting with Leaflet.
>
>
>> *Sharing the app/offline interactivity*
>> I am planning to share this with people who likely do not have R installed
>> on their laptops nor have they ever coded. Does anyone have any
>> suggestions
>> for the best way to do this while retaining interactivity?
>>
>>  Here's the big win of creating a standalone web map. You only have to
> distribute the HTML/CSS/JS and they can be viewed directly (or you also
> supply a tiny server that runs locally and only has to feed the files on a
> localhost port). No need to have a shiny server anywhere, or install R. Its
> simple and clean. It also needs no network connectivity, but you'll not get
> a base map - but you could include a low or medium resolution basemap
> raster in your package.
>
> The only reason to need Shiny here would be if you wanted people to do
> something computational, like click on a bunch of polygons and then fit a
> linear model to the selection, since that would require a round-trip to the
> server for R to compute the fit. (although I suspect there's a JS package
> for linear modelling.... you can do ML in JS these days...)
>
>
>
>> Thank you all, any insight is greatly appreciated.
>>
>> Best,
>> Erin
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>
        [[alternative HTML version deleted]]

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

Re: Adding great circle routes as polylines in Leaflet

Wed, 09/05/2018 - 09:38
Thanks Kent, this looks great!
Hopefully I should be able to convince my client to use this as a
visualization rather than the rather sparse looking current visualization
on a leaflet map.
Regards

Dhiraj Khanna
Mob:09873263331


On Tue, Sep 4, 2018 at 3:55 PM Kent Johnson <[hidden email]> wrote:

>
>
> On Tue, Sep 4, 2018 at 12:18 AM, Dhiraj Khanna <[hidden email]>
> wrote:
>
>> @Kent Johnson <[hidden email]> yes, this is shipping data and you are right, great circle routes are not the best visualization. The problem I am facing is that I have oil trade happening over a period of time from various ports that I need to visualize. Over the selected period of time, there are hundreds of voyages being undertaken by ships. Plotting them all as gc routes looks ugly. My approach has been to classify these ports into regions, which I drew using mapedit and saved them as SF polygons. I then calculated their centroids and those are the coordinates in the ByRoute dataframe. Would appreciate your comments on any other visualization which you think might be appropriate.
>>
>> If your primary interest is to visualize the flows (not the geography)
> then a circle plot might work well. Here is an example showing migration
> flows:
>
> https://gjabel.wordpress.com/2016/05/18/updated-circular-plots-for-directional-bilateral-migration-data/
> made with
> https://cran.r-project.org/web/packages/circlize/index.html
>
> Kent
>
>
>> Regards
>> Dhiraj Khanna
>>
>
        [[alternative HTML version deleted]]

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

Pages