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: 31 min 40 sec ago

Re: correct use of autoKrige and gstat with Digital Elevation Models

Fri, 08/17/2018 - 08:56
Hi Zivan.
You are right. Yes, I did it:

Inverse_Distance_rainfall <- gstat(formula=Cumulata~ExtractedElevationValues, data=df_prec, locations=newgrid)

My problem is when I try to put the output in a raster. How can I do it?
Inverse_Distance_rainfall is a list, but within it I am not able to find the name to use within raster().
Any suggestion?

Thank you
Stefano




         (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----------------
________________________________________
Da: Zivan Karaman [[hidden email]]
Inviato: venerdì 17 agosto 2018 14.50
A: Stefano Sofia
Oggetto: Re: [R-sig-Geo] correct use of autoKrige and gstat with Digital Elevation Models

Have you tried:
Inverse_Distance_rainfall <- gstat(formula=Cumulata~ExtractedElevationValues, data=df_prec, locations=newgrid)
formula, data and locations are NOT the default first 3 arguments to gstat function - see help(gstat).
HTH
Zivan

> Le 16 août 2018 à 09:21, Stefano Sofia <[hidden email]> a écrit :
>
> Dear R-sig-geo list users,
> I am dealing with rainfall isohyet maps, I am using a DEM in order to get more accurate values.
> In particular I am trying to compare Kriging method (through autoKrige) with the Inverse Distance method (through gstat).
> Kriging code works fine, I think that I am not able to make a correct use of gstat.
>
> In both cases,
> - I load my rainfall cumulates in a data frame called df_prec:
>
> Codice_sensore Cumulata Long_Cent Lat_Cent Quota
> 3134 177 13.22888 42.99361  1352
> 3135 98 13.19046 42.85843  1496
> 3136 40 12.81444 43.23840  1005
> 3137 201 13.30681 42.84369  1036
> 3138 145 13.32638 42.90215   980
> 3139 45 13.11180 42.79578  1575
> 3140 69 13.18538 42.91439  1800
> 3234 129 13.44506 42.74417   960
>
> - I load the shape file and the DEM:
>
> ## LOAD SHAPEFILE
> Shape_area5 <- readOGR(dsn="area5.shp", layer="area5")
> Shape_area5_projection <- "+init=epsg:32633 +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
> Shape_area5 <- spTransform(Shape_area5, CRS(Shape_area5_projection))
> ## LOADING DEM
> ita_DEM <- getData('alt', country='ITA', mask=FALSE)
> ita_DEM <- projectRaster(ita_DEM, crs = CRS("+init=epsg:32633"))
>
> - I put the coordinates in the correct way, I evaluate the extracted elevation values and I create a new grid with these new values:
>
> ##
> coordinates(df_prec) <- ~Long_Cent+Lat_Cent
> proj4string(df_prec) <- CRS("+init=epsg:4326")
> df_prec <- spTransform(df_prec, CRS("+init=epsg:32633"))
> ## EXTRACT THE ELEVATION VALUES TO MY POINTS; extract EXPECTS A RASTER INPUT
> df_prec$ExtractedElevationValues <- extract(x=ita_DEM, y=df_prec)
> ## CREATE A NEW GRID
> newgrid <- as(ita_DEM, "SpatialGridDataFrame")
> names(newgrid) <- "ExtractedElevationValues"
>
> - when I apply Kriging, everything works fine
> ## KRIGING
> RegressionKriging_rainfall <- autoKrige(Cumulata~ExtractedElevationValues, df_prec, newgrid)
> prec_rast <- raster(RegressionKriging_rainfall$krige_output)
>
> - if I use gstat in the standard way
>
> Inverse_Distance_rainfall <- gstat(formula=Cumulata~1, data=df_prec, locations=newgrid)
>
> everything goes fine, but this code obviously does not take into consideration the elevation points. When I try to use the inverse distance method like Kriging:
> Inverse_Distance_rainfall <- gstat(formula=Cumulata~ExtractedElevationValues, df_prec, newgrid)
>
> I get the following error:
>
> Error in (function (classes, fdef, mtable)  :
>  unable to find an inherited method for function ‘proj4string’ for signature ‘"NULL"’
>  Calls: my_inversedistanceweighted_DEM -> gstat -> proj4string -> <Anonymous>
>
>
> Where is my mistake?
> Could please somebody help me and give me a hint for a correst use of gstat?
> Thank you for your attention and your help
> Stefano
>
>
>         (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 da Libra ESVA ed  risultato non infetto.
> This message was scanned by Libra ESVA and is believed to be clean.
>
>
>    [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
>  https://urlsand.esvalabs.com/?u=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geo&e=52342f8a&h=f7d4a649&f=y&p=y
--

Questo messaggio  stato analizzato con Libra ESVA ed  risultato non infetto.


________________________________

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 da Libra ESVA ed  risultato non infetto.
This message was scanned 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

Re: st_distance to replace dist in pipes - help needed

Thu, 08/16/2018 - 19:09
Maybe like this?

dist.utm <- matrix(c(st_distance(pts.utm.sf %>%
filter(Attribute=="start")),
                     st_distance(pts.utm.sf %>% filter(Attribute=="end"))),
                     nrow=20,ncol=20)
all.equal(dist.utm, dist.mat)

On Thu, Aug 16, 2018 at 4:53 PM Bannar-Martin, Katherine <
[hidden email]> wrote:

> Hello!
>
> Problem:
> I need to take the following code and replace dist in the "created
> distance matrix" pipe stream (producing object 'res') with st_distance. I'm
> stuck on how to play with the geometry section so that it reproduces the
> same distance matrix as dist.mat. (The numbers will be different because of
> the calculations employed, I know, but the structure of the matrices should
> be the same). Any help is greatly appreciated. Thanks! Katherine
>
> Code using stats::dist:
>
> ###Dependencies####
> require(tidyverse)
> require(sf)
>
> ####Create sample data in decimal degrees:####
>
> START_LONG<-round(runif(n=10,min=-135.7, max=-123.3),digits=4)
> START_LAT<-round(runif(n=10,min=48.26, max=55.80),digits=4)
> END_LONG<-round(runif(n=10,min=-135.7, max=-123.3),digits=4)
> END_LAT<-round(runif(n=10,min=48.28, max=55.82),digits=4)
> df<-as.data.frame(cbind(START_LONG,START_LAT,END_LONG,END_LAT))
>
> ###convert to utm with package sf
> xy1<-df[,c("START_LONG","START_LAT")] #start locations
> colnames(xy1)<-c("X", "Y")
> xy2<-df[,c("END_LONG","END_LAT")] #end locations
> colnames(xy2)<-c("X", "Y")
> xy12 <- rbind(as.matrix(xy1), as.matrix(xy2))
> df12 <- data.frame(data.frame(xy12, Attribute =
> rep(c("start","end"),each=nrow(df)), ID = rep(1:(nrow(df)),2)))
>
> df.SP <- st_as_sf(df12, coords = c("X", "Y"), crs = 4326) #wgs84 long lat
> df.SP<-st_transform(x = df.SP, crs = 32609) #convert to wgs84 utm zone 9
> df.SP$utm_x<-st_coordinates(df.SP)[,1] # get utm coordinates
> df.SP$utm_y<-st_coordinates(df.SP)[,2] # get utm coordinates
> df.SP<-st_set_geometry(df.SP, NULL)# coerce back to data.frame
>
> ###create distance matrix###
> res<- df.SP %>%
>   dplyr::group_by(Attribute) %>%
>   mutate(rows = row_number()) %>%
>   left_join(df.SP, by = c('Attribute')) %>%
>   rowwise() %>%
>   mutate(Distance = ifelse(dist(rbind(c(utm_x.x, utm_y.x), c(utm_x.y,
> utm_y.y))) != 0,
>                            dist(rbind(c(utm_x.x, utm_y.x), c(utm_x.y,
> utm_y.y))), NA))
>
> dist.mat<-matrix(res$Distance,nrow = 20,ncol = 20)
>
> Data prep for eventually using sf::st_distance:
> pts.utm.sf <- df12 %>%
>   st_as_sf(coords=c('X','Y'), crs=4326) %>% #wgs84 long lat
>   st_transform(32609) #wgs84 utm zone9
>
> #distance matrix without groups (note using group_by "Attribute" does not
> get me the matrix I want):
> dist.utm<- pts.utm.sf %>% st_distance() #calculate distances on utm in m
>
>         [[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

st_distance to replace dist in pipes - help needed

Thu, 08/16/2018 - 15:52
Hello!

Problem:
I need to take the following code and replace dist in the "created distance matrix" pipe stream (producing object 'res') with st_distance. I'm stuck on how to play with the geometry section so that it reproduces the same distance matrix as dist.mat. (The numbers will be different because of the calculations employed, I know, but the structure of the matrices should be the same). Any help is greatly appreciated. Thanks! Katherine

Code using stats::dist:

###Dependencies####
require(tidyverse)
require(sf)

####Create sample data in decimal degrees:####

START_LONG<-round(runif(n=10,min=-135.7, max=-123.3),digits=4)
START_LAT<-round(runif(n=10,min=48.26, max=55.80),digits=4)
END_LONG<-round(runif(n=10,min=-135.7, max=-123.3),digits=4)
END_LAT<-round(runif(n=10,min=48.28, max=55.82),digits=4)
df<-as.data.frame(cbind(START_LONG,START_LAT,END_LONG,END_LAT))

###convert to utm with package sf
xy1<-df[,c("START_LONG","START_LAT")] #start locations
colnames(xy1)<-c("X", "Y")
xy2<-df[,c("END_LONG","END_LAT")] #end locations
colnames(xy2)<-c("X", "Y")
xy12 <- rbind(as.matrix(xy1), as.matrix(xy2))
df12 <- data.frame(data.frame(xy12, Attribute = rep(c("start","end"),each=nrow(df)), ID = rep(1:(nrow(df)),2)))

df.SP <- st_as_sf(df12, coords = c("X", "Y"), crs = 4326) #wgs84 long lat
df.SP<-st_transform(x = df.SP, crs = 32609) #convert to wgs84 utm zone 9
df.SP$utm_x<-st_coordinates(df.SP)[,1] # get utm coordinates
df.SP$utm_y<-st_coordinates(df.SP)[,2] # get utm coordinates
df.SP<-st_set_geometry(df.SP, NULL)# coerce back to data.frame

###create distance matrix###
res<- df.SP %>%
  dplyr::group_by(Attribute) %>%
  mutate(rows = row_number()) %>%
  left_join(df.SP, by = c('Attribute')) %>%
  rowwise() %>%
  mutate(Distance = ifelse(dist(rbind(c(utm_x.x, utm_y.x), c(utm_x.y, utm_y.y))) != 0,
                           dist(rbind(c(utm_x.x, utm_y.x), c(utm_x.y, utm_y.y))), NA))

dist.mat<-matrix(res$Distance,nrow = 20,ncol = 20)

Data prep for eventually using sf::st_distance:
pts.utm.sf <- df12 %>%
  st_as_sf(coords=c('X','Y'), crs=4326) %>% #wgs84 long lat
  st_transform(32609) #wgs84 utm zone9

#distance matrix without groups (note using group_by "Attribute" does not get me the matrix I want):
dist.utm<- pts.utm.sf %>% st_distance() #calculate distances on utm in m

        [[alternative HTML version deleted]]

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

Morans I

Thu, 08/16/2018 - 06:44
hi all,

is the below the correct way to access spatial autocorrelation using morans I in a glm:

#get coordinates first
coords<-as.matrix(cbind(data$long,data$lat))

#model
m1 <- glm(response~ predictor1 + predictor 2+ predictor1*predictor2, family=binomial, data=data)

# Compute Moran's I using residuals of model
lstw <- nb2listw((knn2nb(knearneigh(coords, k=1))))

moran.test(residuals(m1), lstw)


Advise from experts using function moran.test will be much appreciated.
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

correct use of autoKrige and gstat with Digital Elevation Models

Thu, 08/16/2018 - 02:21
Dear R-sig-geo list users,
I am dealing with rainfall isohyet maps, I am using a DEM in order to get more accurate values.
In particular I am trying to compare Kriging method (through autoKrige) with the Inverse Distance method (through gstat).
Kriging code works fine, I think that I am not able to make a correct use of gstat.

In both cases,
- I load my rainfall cumulates in a data frame called df_prec:

Codice_sensore Cumulata Long_Cent Lat_Cent Quota
3134 177 13.22888 42.99361  1352
3135 98 13.19046 42.85843  1496
3136 40 12.81444 43.23840  1005
3137 201 13.30681 42.84369  1036
3138 145 13.32638 42.90215   980
3139 45 13.11180 42.79578  1575
3140 69 13.18538 42.91439  1800
3234 129 13.44506 42.74417   960

- I load the shape file and the DEM:

## LOAD SHAPEFILE
Shape_area5 <- readOGR(dsn="area5.shp", layer="area5")
Shape_area5_projection <- "+init=epsg:32633 +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
Shape_area5 <- spTransform(Shape_area5, CRS(Shape_area5_projection))
## LOADING DEM
ita_DEM <- getData('alt', country='ITA', mask=FALSE)
ita_DEM <- projectRaster(ita_DEM, crs = CRS("+init=epsg:32633"))

- I put the coordinates in the correct way, I evaluate the extracted elevation values and I create a new grid with these new values:

##
coordinates(df_prec) <- ~Long_Cent+Lat_Cent
proj4string(df_prec) <- CRS("+init=epsg:4326")
df_prec <- spTransform(df_prec, CRS("+init=epsg:32633"))
## EXTRACT THE ELEVATION VALUES TO MY POINTS; extract EXPECTS A RASTER INPUT
df_prec$ExtractedElevationValues <- extract(x=ita_DEM, y=df_prec)
## CREATE A NEW GRID
newgrid <- as(ita_DEM, "SpatialGridDataFrame")
names(newgrid) <- "ExtractedElevationValues"

- when I apply Kriging, everything works fine
## KRIGING
RegressionKriging_rainfall <- autoKrige(Cumulata~ExtractedElevationValues, df_prec, newgrid)
prec_rast <- raster(RegressionKriging_rainfall$krige_output)

- if I use gstat in the standard way

Inverse_Distance_rainfall <- gstat(formula=Cumulata~1, data=df_prec, locations=newgrid)

everything goes fine, but this code obviously does not take into consideration the elevation points. When I try to use the inverse distance method like Kriging:
Inverse_Distance_rainfall <- gstat(formula=Cumulata~ExtractedElevationValues, df_prec, newgrid)

I get the following error:

Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function �proj4string� for signature �"NULL"�
  Calls: my_inversedistanceweighted_DEM -> gstat -> proj4string -> <Anonymous>


Where is my mistake?
Could please somebody help me and give me a hint for a correst use of gstat?
Thank you for your attention and your help
Stefano


         (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 da Libra ESVA ed  risultato non infetto.
This message was scanned by Libra ESVA and is believed to be clean.


        [[alternative HTML version deleted]]


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

Re: how to create several polygons from a list of vertices

Wed, 08/15/2018 - 18:36
Hi Vijay, Michael and Obrl

Thanks for the answers, finally I could do the job.

I learned a lot too!

Best regards

Antonio Olinto

Em Ter, 14 de ago de 2018 11:19 PM, Vijay Lulla <[hidden email]>
escreveu:

> Maybe you can try this then?
>
> polys <- SpatialPolygons(lapply(unique(vertices$cod),
>                          function(x) {
>                            Polygons(list(Polygon(vertices[vertices$cod ==
> x, 1:2])), ID=x)
>                          } ))
> HTH,
> Vijay.
>
> On Tue, Aug 14, 2018 at 6:03 PM Antonio Silva <[hidden email]>
> wrote:
>
>> Thanks Lulla,
>>
>> Nice solution. I could also export it as a shapefile after transforming
>> it to a spatial polygon dataframe.
>>
>> The problem is that I could not "individualize" the squares in a
>> multipart layer. They all have the same ID. I tried to change this without
>> success: "Single ID required".
>>
>> The attribute table of the shapefile should have 6 lines in my example
>> and not only one.
>>
>> Any other option?
>>
>> Thanks again,
>>
>> Antonio Olinto
>>
>>
>> Em ter, 14 de ago de 2018 às 18:10, Vijay Lulla <[hidden email]>
>> escreveu:
>>
>>> Maybe something like this?
>>>
>>> poly <- SpatialPolygons(list(Polygons(tapply(seq_len(nrow(vertices)),
>>>                                              vertices$cod,
>>>                                              function(x)
>>> Polygon(vertices[x,1:2])), ID="1")),
>>>                         proj4string=CRS("+proj=longlat +ellps=WGS84
>>> +datum=WGS84 +no_defs"))
>>>
>>>
>>> On Tue, Aug 14, 2018 at 4:17 PM Antonio Silva <[hidden email]>
>>> wrote:
>>>
>>>> Hi,
>>>>
>>>> I have a data.frame with the vertices (lon / lat) and codes from several
>>>> squares (more than 500 in the real dataset).
>>>> I want to create an object with these polygons (squares) and after this
>>>> export it as a shapefile.
>>>> With the script below I can draw one square.
>>>> library(sp)
>>>> P1 = Polygon(vertices[1:4,1:2])
>>>> Ps1 = SpatialPolygons(list(Polygons(list(P1), ID = "1")),
>>>> proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
>>>> plot(Ps1, axes = TRUE)
>>>>
>>>> Now I'm trying to create one object with all squares at once.
>>>> Is it possible?
>>>>
>>>> Thanks a lot,
>>>>
>>>> Antônio Olinto
>>>>
>>>> sample data:vertices
>>>>    lon lat cod
>>>> 1  -33 -23   1
>>>> 2  -32 -23   1
>>>> 3  -32 -22   1
>>>> 4  -33 -22   1
>>>> 5  -32 -23   2
>>>> 6  -31 -23   2
>>>> 7  -31 -22   2
>>>> 8  -32 -22   2
>>>> 9  -31 -23   3
>>>> 10 -30 -23   3
>>>> 11 -30 -22   3
>>>> 12 -31 -22   3
>>>> 13 -33 -22   4
>>>> 14 -32 -22   4
>>>> 15 -32 -21   4
>>>> 16 -33 -21   4
>>>> 17 -32 -22   5
>>>> 18 -31 -22   5
>>>> 19 -31 -21   5
>>>> 20 -32 -21   5
>>>> 21 -31 -22   6
>>>> 22 -30 -22   6
>>>> 23 -30 -21   6
>>>> 24 -31 -21   6
>>>>
>>>>         [[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: [localmoran.exact]

Wed, 08/15/2018 - 05:50

> Le 14 août 2018 à 20:09, Roger Bivand <[hidden email]> a écrit :
>
> If you post, you then reply in-thread on the list. Otherwise, nobody else knows what the resolution was. :
Sorry. My bad.
>
> Your explanation below contains lots of inaccuracies - all Spatial*DataFrames have standard data.frame access methods as I showed.
No doubt about it ! I will try do better now that I know about it.  

Thank you Sir.
Lisa

>
> Roger Bivand
>
> On Tue, 14 Aug 2018, Lisa Menez wrote:
>
>> Sir Bivand,
>>
>> The complete example was here :
>>
>> My database is available here : http://hcc.cnrs.fr/wp-content/uploads/2018/08/mydata2018-1.txt <http://hcc.cnrs.fr/wp-content/uploads/2018/08/mydata2018-1.txt> and my code : http://hcc.cnrs.fr/wp-content/uploads/2018/08/myCode_2018.txt <http://hcc.cnrs.fr/wp-content/uploads/2018/08/myCode_2018.txt>
>>
>> The localmoran.sad() actually works. Is it because « If we assume global spatial independence, we can give for any projection matrix M the spectrum of eigenvalues as well as the saddle point approximation of local Moran’s Ii in analytical terms. This increases the computational efficiency as it avoids reverting to numerical calculations of the spectrum of eigenvalues and finding iteratively the root of the saddle point equation » ?
>>
>> To the best of my knowledge, I am constrained by the settings of my data (extracted from Eurostat) on calling EU_NUTS@data rather than EU_NUTS. Indeed EU_NUTS are spatial polygons and EU_NUTS@data their relative attributes. I failed to change this name.
>>
>> for each i, EU_NUTS@data[,year[i]] is a numerical vector. year is a character vector containing the sample years, allowing to go through the corresponding EU_NUTS@data' columns.
>>
>> On your advice, I tried to go through debug(localmoran.exact). As I could not find any good reason why (2*t2<t1^2), I turned to your last remark : Changing NUTS_L1 into a factor before using it into the loop. It actually was the issue. I miscoded the factor, keeping character chains and not transforming them into numbers.
>>
>>
>> Thank you very much for your help and recommandations.
>> Lisa
>>
>>
>>> Le 13 août 2018 à 16:12, Roger Bivand <[hidden email]> a écrit :
>>>
>>> On Mon, 13 Aug 2018, Lisa Menez wrote:
>>>
>>>>
>>>> Dear all,
>>>>
>>>> I have been trying to obtain local Moran’s I with the exact approach (localmoran.exact) but I have some issues and would need your help.
>>>
>>> A complete reproducible example is required (script and data, data may be downloaded, do not attach). Does the same thing happen with localmoran.sad()?
>>>
>>>>
>>>> I am trying to bring to light how spatial patterns of wealth have evolved from 2000 to 2016 in the European Union.
>>>>
>>>> Looking at the Regional Gross Domestic Product (GDP) of European regions, my lm model looks like this :
>>>>
>>>> lm(log(EU_NUTS@data[,year[i]])~1+factor(EU_NUTS@data$NUTS_L1))
>>>>
>>>> where log(EU_NUTS@data[,year[i]]) are log of regional GDPs
>>>>
>>>
>>> Never, ever use @ to access data. EU_NUTS@data$NUTS_L1 should be EU_NUTS$NUTS_L1 and should be (made a) factor first. It isn't clear what log(EU_NUTS@data[,year[i]]) will end up being - is "year" a character vector? If so, EU_NUTS[[year[i]]] is the correct syntax.
>>>
>>>> and factor(EU_NUTS@data$NUTS_L1) is the vector of country dummies.
>>>>
>>>> Sadly, I obtain an error message that says :
>>>>
>>>> Error in integrate(integrand, lower = 0, upper = upper) :
>>>> a limit is missing
>>>> Moreover : Warning messages:
>>>> 1: In sqrt(2 * t2 - t1^2) : production of NaN
>>>> 2: In sqrt(2 * t2 - t1^2) : production of NaN
>>>>
>>>> According to the GitHub source code ,
>>>
>>> No, as you can see from https://cran.r-project.org/package=spdep, the development site of spdep is https://github.com/r-spatial/spdep/. Never use non-authorised copies. To see the source, simply type the function name. Use debug(localmoran.exact) to step through the function while it runs, to see the values of the objects you are in doubt about. Do not do this in RStudio, its debugging interface tries to be far too clever.
>>>
>>> The following is illegible in plain text; post in plain text only.
>>>
>>> Roger
>>>
>>>>
>>>>      Vi <- listw2star <https://rdrr.io/cran/spdep/man/localmoran.sad.html>(B, select[i], style=style, n, D <https://rdrr.io/r/stats/deriv.html>, a,
>>>>    zero.policy=zero.policy)
>>>>      Viu <- lag.listw <https://rdrr.io/cran/spdep/man/lag.listw.html>(Vi, u, zero.policy=TRUE <https://rdrr.io/r/base/logical.html>)
>>>> Ii <- c <https://rdrr.io/r/base/c.html>(crossprod <https://rdrr.io/r/base/crossprod.html>(u, Viu) / utu)
>>>>      ViX <- lag.listw <https://rdrr.io/cran/spdep/man/lag.listw.html>(Vi, X, zero.policy=TRUE <https://rdrr.io/r/base/logical.html>)
>>>>      MViM <- t <https://rdrr.io/r/base/t.html>(X) %*% ViX %*% XtXinv
>>>>      t1 <- -sum <https://rdrr.io/r/base/sum.html>(diag <https://rdrr.io/r/base/diag.html>(MViM))
>>>>      sumsq.Vi <- function <https://rdrr.io/r/base/function.html>(x) {
>>>>          if <https://rdrr.io/r/base/Control.html> (is.null <https://rdrr.io/r/base/NULL.html>(x)) NA <https://rdrr.io/r/base/NA.html>
>>>>    else <https://rdrr.io/r/base/Control.html> sum <https://rdrr.io/r/base/sum.html>(x^2)
>>>> }
>>>> trVi2 <- sum <https://rdrr.io/r/base/sum.html>(sapply <https://rdrr.io/r/base/lapply.html>(Vi$weights <https://rdrr.io/r/stats/weights.html>, sumsq.Vi), na.rm=TRUE <https://rdrr.io/r/base/logical.html>)
>>>> t2a <- sum <https://rdrr.io/r/base/sum.html>(diag <https://rdrr.io/r/base/diag.html>(t <https://rdrr.io/r/base/t.html>(ViX) %*% ViX %*% XtXinv))
>>>> t2b <- sum <https://rdrr.io/r/base/sum.html>(diag <https://rdrr.io/r/base/diag.html>(MViM %*% MViM))
>>>> t2 <- trVi2 - 2*t2a + t2b
>>>> e1 <- 0.5 * (t1 + sqrt <https://rdrr.io/r/base/MathFun.html>(2*t2 - t1^2))
>>>> en <- 0.5 * (t1 - sqrt <https://rdrr.io/r/base/MathFun.html>(2*t2 - t1^2))
>>>>
>>>> My guess is that my setting causes (2*t2<t1^2) : Is it ? Why would it ?
>>>>
>>>> My database is available here : http://hcc.cnrs.fr/wp-content/uploads/2018/08/mydata2018-1.txt <http://hcc.cnrs.fr/wp-content/uploads/2018/08/mydata2018-1.txt>
>>>> and my code : http://hcc.cnrs.fr/wp-content/uploads/2018/08/myCode_2018.txt <http://hcc.cnrs.fr/wp-content/uploads/2018/08/myCode_2018.txt>
>>>>
>>>> I thank you very much for your attention and your help.
>>>>
>>>> Best Regards
>>>> Lisa
>>>>
>>>> PhD Student
>>>> GREDEG, Université Côté d’Azur
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> [[alternative HTML version deleted]]
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> [hidden email]
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>
>>> --
>>> Roger Bivand
>>> Department of Economics, Norwegian School of Economics,
>>> Helleveien 30, N-5045 Bergen, Norway.
>>> voice: +47 55 95 93 55; e-mail: [hidden email]
>>> http://orcid.org/0000-0003-2392-6140
>>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>>
>>
>
> --
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; e-mail: [hidden email]
> http://orcid.org/0000-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: PhD course ECS530, Bergen, 3-8 December

Wed, 08/15/2018 - 04:18
On Wed, 15 Aug 2018, Roger Bivand wrote:

> A PhD-level course in spatial data analysis (ECS530, Auditorium C) will be
> held 3-8 December 2018 in Bergen:
>
> https://www.nhh.no/en/courses/analysing-spatial-data/
>
> External participants should apply using this form:
>
> https://www.nhh.no/en/study-programmes/phd-programme-at-nhh/phd-courses/application-to-follow-phd-courses-at-nhh/
>
> For further details see the course page.
I've already been asked about support: participants must cover their own
travel and living costs; no support is offered, I'm afraid, apart from
free tuition.

Roger

>
>

--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway

PhD course ECS530, Bergen, 3-8 December

Wed, 08/15/2018 - 03:51
A PhD-level course in spatial data analysis (ECS530, Auditorium C) will be
held 3-8 December 2018 in Bergen:

https://www.nhh.no/en/courses/analysing-spatial-data/

External participants should apply using this form:

https://www.nhh.no/en/study-programmes/phd-programme-at-nhh/phd-courses/application-to-follow-phd-courses-at-nhh/

For further details see the course page.

--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway

Re: how to create several polygons from a list of vertices

Tue, 08/14/2018 - 21:19
Maybe you can try this then?

polys <- SpatialPolygons(lapply(unique(vertices$cod),
                         function(x) {
                           Polygons(list(Polygon(vertices[vertices$cod ==
x, 1:2])), ID=x)
                         } ))
HTH,
Vijay.

On Tue, Aug 14, 2018 at 6:03 PM Antonio Silva <[hidden email]> wrote:

> Thanks Lulla,
>
> Nice solution. I could also export it as a shapefile after transforming it
> to a spatial polygon dataframe.
>
> The problem is that I could not "individualize" the squares in a multipart
> layer. They all have the same ID. I tried to change this without success:
> "Single ID required".
>
> The attribute table of the shapefile should have 6 lines in my example and
> not only one.
>
> Any other option?
>
> Thanks again,
>
> Antonio Olinto
>
>
> Em ter, 14 de ago de 2018 às 18:10, Vijay Lulla <[hidden email]>
> escreveu:
>
>> Maybe something like this?
>>
>> poly <- SpatialPolygons(list(Polygons(tapply(seq_len(nrow(vertices)),
>>                                              vertices$cod,
>>                                              function(x)
>> Polygon(vertices[x,1:2])), ID="1")),
>>                         proj4string=CRS("+proj=longlat +ellps=WGS84
>> +datum=WGS84 +no_defs"))
>>
>>
>> On Tue, Aug 14, 2018 at 4:17 PM Antonio Silva <[hidden email]>
>> wrote:
>>
>>> Hi,
>>>
>>> I have a data.frame with the vertices (lon / lat) and codes from several
>>> squares (more than 500 in the real dataset).
>>> I want to create an object with these polygons (squares) and after this
>>> export it as a shapefile.
>>> With the script below I can draw one square.
>>> library(sp)
>>> P1 = Polygon(vertices[1:4,1:2])
>>> Ps1 = SpatialPolygons(list(Polygons(list(P1), ID = "1")),
>>> proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
>>> plot(Ps1, axes = TRUE)
>>>
>>> Now I'm trying to create one object with all squares at once.
>>> Is it possible?
>>>
>>> Thanks a lot,
>>>
>>> Antônio Olinto
>>>
>>> sample data:vertices
>>>    lon lat cod
>>> 1  -33 -23   1
>>> 2  -32 -23   1
>>> 3  -32 -22   1
>>> 4  -33 -22   1
>>> 5  -32 -23   2
>>> 6  -31 -23   2
>>> 7  -31 -22   2
>>> 8  -32 -22   2
>>> 9  -31 -23   3
>>> 10 -30 -23   3
>>> 11 -30 -22   3
>>> 12 -31 -22   3
>>> 13 -33 -22   4
>>> 14 -32 -22   4
>>> 15 -32 -21   4
>>> 16 -33 -21   4
>>> 17 -32 -22   5
>>> 18 -31 -22   5
>>> 19 -31 -21   5
>>> 20 -32 -21   5
>>> 21 -31 -22   6
>>> 22 -30 -22   6
>>> 23 -30 -21   6
>>> 24 -31 -21   6
>>>
>>>         [[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: how to create several polygons from a list of vertices

Tue, 08/14/2018 - 19:09
Here's another way with spbabel.


library(dplyr)
library(spbabel)
pts <-
  tibble::tribble(~ID,  ~x,  ~y, ~grp,
                   1 , -33, -23,    1,
                   2 , -32, -23,    1,
                   3 , -32, -22,    1,
                   4 , -33, -22,    1,
                   5 , -32, -23,    2,
                   6 , -31, -23,    2,
                   7 , -31, -22,    2,
                   8 , -32, -22,    2,
                   9 , -31, -23,    3,
                   10, -30, -23,    3,
                   11, -30, -22,    3,
                   12, -31, -22,    3,
                   13, -33, -22,    4,
                   14, -32, -22,    4,
                   15, -32, -21,    4,
                   16, -33, -21,    4,
                   17, -32, -22,    5,
                   18, -31, -22,    5,
                   19, -31, -21,    5,
                   20, -32, -21,    5,
                   21, -31, -22,    6,
                   22, -30, -22,    6,
                   23, -30, -21,    6,
                   24, -31, -21,    6)


## objects and branches (parts) are the same level

## objects and branches (parts) are the same level
## and we maintain "grp" as the object attribute data
x <- pts %>% transmute(grp = grp, x_ = x, y_ = y,
               object_ = grp, branch_ = grp,
               order_ = ID, island_ = TRUE) %>% spbabel::sp()





On Wed, 15 Aug 2018 at 09:25 obrl soil <[hidden email]> wrote:

> Hi Antonio,
>
> have you tried with sf? Like:
>
> library(sf)
>
> pts <-
>   tibble::tribble(~ID,  ~x,  ~y, ~grp,
>                    1 , -33, -23,    1,
>                    2 , -32, -23,    1,
>                    3 , -32, -22,    1,
>                    4 , -33, -22,    1,
>                    5 , -32, -23,    2,
>                    6 , -31, -23,    2,
>                    7 , -31, -22,    2,
>                    8 , -32, -22,    2,
>                    9 , -31, -23,    3,
>                    10, -30, -23,    3,
>                    11, -30, -22,    3,
>                    12, -31, -22,    3,
>                    13, -33, -22,    4,
>                    14, -32, -22,    4,
>                    15, -32, -21,    4,
>                    16, -33, -21,    4,
>                    17, -32, -22,    5,
>                    18, -31, -22,    5,
>                    19, -31, -21,    5,
>                    20, -32, -21,    5,
>                    21, -31, -22,    6,
>                    22, -30, -22,    6,
>                    23, -30, -21,    6,
>                    24, -31, -21,    6)
>
> squares <- split(pts, pts$grp)
> squares <- lapply(squares, function(g) {
>   # get just coords
>   g <- as.matrix(g[, c(2,3)])
>   # repeat first point last to close poly
>   g <- rbind(g, g[1, ])
>   # convert to an sf polygon object
>   gp <- sf::st_polygon(list(g))
>   # make sure the vertices are in an order
>   # that complies with the simple
>   # features standard
>   gp <- sf::st_buffer(gp, 0L)
> })
> # turn list of polygons into geometry column
> squares_sfc <- sf::st_sfc(squares) # can add crs = ?? to this call
> # add an ID to make an sf data frame
> squares_sf <- sf::st_sf('ID' = seq(6), 'geometry' = squares_sfc)
>
> # if you still need to use sp for whatever reason
> squares_sp <- as(squares_sf, "Spatial")
>
> Regards,
> @obrl_soil
>
> On Wed, Aug 15, 2018 at 8:03 AM, Antonio Silva <[hidden email]>
> wrote:
> > Thanks Lulla,
> >
> > Nice solution. I could also export it as a shapefile after transforming
> it
> > to a spatial polygon dataframe.
> >
> > The problem is that I could not "individualize" the squares in a
> multipart
> > layer. They all have the same ID. I tried to change this without success:
> > "Single ID required".
> >
> > The attribute table of the shapefile should have 6 lines in my example
> and
> > not only one.
> >
> > Any other option?
> >
> > Thanks again,
> >
> > Antonio Olinto
> >
> >
> > Em ter, 14 de ago de 2018 às 18:10, Vijay Lulla <[hidden email]>
> > escreveu:
> >
> >> Maybe something like this?
> >>
> >> poly <- SpatialPolygons(list(Polygons(tapply(seq_len(nrow(vertices)),
> >>                                              vertices$cod,
> >>                                              function(x)
> >> Polygon(vertices[x,1:2])), ID="1")),
> >>                         proj4string=CRS("+proj=longlat +ellps=WGS84
> >> +datum=WGS84 +no_defs"))
> >>
> >>
> >> On Tue, Aug 14, 2018 at 4:17 PM Antonio Silva <[hidden email]>
> >> wrote:
> >>
> >>> Hi,
> >>>
> >>> I have a data.frame with the vertices (lon / lat) and codes from
> several
> >>> squares (more than 500 in the real dataset).
> >>> I want to create an object with these polygons (squares) and after this
> >>> export it as a shapefile.
> >>> With the script below I can draw one square.
> >>> library(sp)
> >>> P1 = Polygon(vertices[1:4,1:2])
> >>> Ps1 = SpatialPolygons(list(Polygons(list(P1), ID = "1")),
> >>> proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
> >>> plot(Ps1, axes = TRUE)
> >>>
> >>> Now I'm trying to create one object with all squares at once.
> >>> Is it possible?
> >>>
> >>> Thanks a lot,
> >>>
> >>> Antônio Olinto
> >>>
> >>> sample data:vertices
> >>>    lon lat cod
> >>> 1  -33 -23   1
> >>> 2  -32 -23   1
> >>> 3  -32 -22   1
> >>> 4  -33 -22   1
> >>> 5  -32 -23   2
> >>> 6  -31 -23   2
> >>> 7  -31 -22   2
> >>> 8  -32 -22   2
> >>> 9  -31 -23   3
> >>> 10 -30 -23   3
> >>> 11 -30 -22   3
> >>> 12 -31 -22   3
> >>> 13 -33 -22   4
> >>> 14 -32 -22   4
> >>> 15 -32 -21   4
> >>> 16 -33 -21   4
> >>> 17 -32 -22   5
> >>> 18 -31 -22   5
> >>> 19 -31 -21   5
> >>> 20 -32 -21   5
> >>> 21 -31 -22   6
> >>> 22 -30 -22   6
> >>> 23 -30 -21   6
> >>> 24 -31 -21   6
> >>>
> >>>         [[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
> --
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: how to create several polygons from a list of vertices

Tue, 08/14/2018 - 18:25
Hi Antonio,

have you tried with sf? Like:

library(sf)

pts <-
  tibble::tribble(~ID,  ~x,  ~y, ~grp,
                   1 , -33, -23,    1,
                   2 , -32, -23,    1,
                   3 , -32, -22,    1,
                   4 , -33, -22,    1,
                   5 , -32, -23,    2,
                   6 , -31, -23,    2,
                   7 , -31, -22,    2,
                   8 , -32, -22,    2,
                   9 , -31, -23,    3,
                   10, -30, -23,    3,
                   11, -30, -22,    3,
                   12, -31, -22,    3,
                   13, -33, -22,    4,
                   14, -32, -22,    4,
                   15, -32, -21,    4,
                   16, -33, -21,    4,
                   17, -32, -22,    5,
                   18, -31, -22,    5,
                   19, -31, -21,    5,
                   20, -32, -21,    5,
                   21, -31, -22,    6,
                   22, -30, -22,    6,
                   23, -30, -21,    6,
                   24, -31, -21,    6)

squares <- split(pts, pts$grp)
squares <- lapply(squares, function(g) {
  # get just coords
  g <- as.matrix(g[, c(2,3)])
  # repeat first point last to close poly
  g <- rbind(g, g[1, ])
  # convert to an sf polygon object
  gp <- sf::st_polygon(list(g))
  # make sure the vertices are in an order
  # that complies with the simple
  # features standard
  gp <- sf::st_buffer(gp, 0L)
})
# turn list of polygons into geometry column
squares_sfc <- sf::st_sfc(squares) # can add crs = ?? to this call
# add an ID to make an sf data frame
squares_sf <- sf::st_sf('ID' = seq(6), 'geometry' = squares_sfc)

# if you still need to use sp for whatever reason
squares_sp <- as(squares_sf, "Spatial")

Regards,
@obrl_soil

On Wed, Aug 15, 2018 at 8:03 AM, Antonio Silva <[hidden email]> wrote:
> Thanks Lulla,
>
> Nice solution. I could also export it as a shapefile after transforming it
> to a spatial polygon dataframe.
>
> The problem is that I could not "individualize" the squares in a multipart
> layer. They all have the same ID. I tried to change this without success:
> "Single ID required".
>
> The attribute table of the shapefile should have 6 lines in my example and
> not only one.
>
> Any other option?
>
> Thanks again,
>
> Antonio Olinto
>
>
> Em ter, 14 de ago de 2018 às 18:10, Vijay Lulla <[hidden email]>
> escreveu:
>
>> Maybe something like this?
>>
>> poly <- SpatialPolygons(list(Polygons(tapply(seq_len(nrow(vertices)),
>>                                              vertices$cod,
>>                                              function(x)
>> Polygon(vertices[x,1:2])), ID="1")),
>>                         proj4string=CRS("+proj=longlat +ellps=WGS84
>> +datum=WGS84 +no_defs"))
>>
>>
>> On Tue, Aug 14, 2018 at 4:17 PM Antonio Silva <[hidden email]>
>> wrote:
>>
>>> Hi,
>>>
>>> I have a data.frame with the vertices (lon / lat) and codes from several
>>> squares (more than 500 in the real dataset).
>>> I want to create an object with these polygons (squares) and after this
>>> export it as a shapefile.
>>> With the script below I can draw one square.
>>> library(sp)
>>> P1 = Polygon(vertices[1:4,1:2])
>>> Ps1 = SpatialPolygons(list(Polygons(list(P1), ID = "1")),
>>> proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
>>> plot(Ps1, axes = TRUE)
>>>
>>> Now I'm trying to create one object with all squares at once.
>>> Is it possible?
>>>
>>> Thanks a lot,
>>>
>>> Antônio Olinto
>>>
>>> sample data:vertices
>>>    lon lat cod
>>> 1  -33 -23   1
>>> 2  -32 -23   1
>>> 3  -32 -22   1
>>> 4  -33 -22   1
>>> 5  -32 -23   2
>>> 6  -31 -23   2
>>> 7  -31 -22   2
>>> 8  -32 -22   2
>>> 9  -31 -23   3
>>> 10 -30 -23   3
>>> 11 -30 -22   3
>>> 12 -31 -22   3
>>> 13 -33 -22   4
>>> 14 -32 -22   4
>>> 15 -32 -21   4
>>> 16 -33 -21   4
>>> 17 -32 -22   5
>>> 18 -31 -22   5
>>> 19 -31 -21   5
>>> 20 -32 -21   5
>>> 21 -31 -22   6
>>> 22 -30 -22   6
>>> 23 -30 -21   6
>>> 24 -31 -21   6
>>>
>>>         [[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: how to create several polygons from a list of vertices

Tue, 08/14/2018 - 17:03
Thanks Lulla,

Nice solution. I could also export it as a shapefile after transforming it
to a spatial polygon dataframe.

The problem is that I could not "individualize" the squares in a multipart
layer. They all have the same ID. I tried to change this without success:
"Single ID required".

The attribute table of the shapefile should have 6 lines in my example and
not only one.

Any other option?

Thanks again,

Antonio Olinto


Em ter, 14 de ago de 2018 às 18:10, Vijay Lulla <[hidden email]>
escreveu:

> Maybe something like this?
>
> poly <- SpatialPolygons(list(Polygons(tapply(seq_len(nrow(vertices)),
>                                              vertices$cod,
>                                              function(x)
> Polygon(vertices[x,1:2])), ID="1")),
>                         proj4string=CRS("+proj=longlat +ellps=WGS84
> +datum=WGS84 +no_defs"))
>
>
> On Tue, Aug 14, 2018 at 4:17 PM Antonio Silva <[hidden email]>
> wrote:
>
>> Hi,
>>
>> I have a data.frame with the vertices (lon / lat) and codes from several
>> squares (more than 500 in the real dataset).
>> I want to create an object with these polygons (squares) and after this
>> export it as a shapefile.
>> With the script below I can draw one square.
>> library(sp)
>> P1 = Polygon(vertices[1:4,1:2])
>> Ps1 = SpatialPolygons(list(Polygons(list(P1), ID = "1")),
>> proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
>> plot(Ps1, axes = TRUE)
>>
>> Now I'm trying to create one object with all squares at once.
>> Is it possible?
>>
>> Thanks a lot,
>>
>> Antônio Olinto
>>
>> sample data:vertices
>>    lon lat cod
>> 1  -33 -23   1
>> 2  -32 -23   1
>> 3  -32 -22   1
>> 4  -33 -22   1
>> 5  -32 -23   2
>> 6  -31 -23   2
>> 7  -31 -22   2
>> 8  -32 -22   2
>> 9  -31 -23   3
>> 10 -30 -23   3
>> 11 -30 -22   3
>> 12 -31 -22   3
>> 13 -33 -22   4
>> 14 -32 -22   4
>> 15 -32 -21   4
>> 16 -33 -21   4
>> 17 -32 -22   5
>> 18 -31 -22   5
>> 19 -31 -21   5
>> 20 -32 -21   5
>> 21 -31 -22   6
>> 22 -30 -22   6
>> 23 -30 -21   6
>> 24 -31 -21   6
>>
>>         [[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: how to create several polygons from a list of vertices

Tue, 08/14/2018 - 16:10
Maybe something like this?

poly <- SpatialPolygons(list(Polygons(tapply(seq_len(nrow(vertices)),
                                             vertices$cod,
                                             function(x)
Polygon(vertices[x,1:2])), ID="1")),
                        proj4string=CRS("+proj=longlat +ellps=WGS84
+datum=WGS84 +no_defs"))


On Tue, Aug 14, 2018 at 4:17 PM Antonio Silva <[hidden email]> wrote:

> Hi,
>
> I have a data.frame with the vertices (lon / lat) and codes from several
> squares (more than 500 in the real dataset).
> I want to create an object with these polygons (squares) and after this
> export it as a shapefile.
> With the script below I can draw one square.
> library(sp)
> P1 = Polygon(vertices[1:4,1:2])
> Ps1 = SpatialPolygons(list(Polygons(list(P1), ID = "1")),
> proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
> plot(Ps1, axes = TRUE)
>
> Now I'm trying to create one object with all squares at once.
> Is it possible?
>
> Thanks a lot,
>
> Antônio Olinto
>
> sample data:vertices
>    lon lat cod
> 1  -33 -23   1
> 2  -32 -23   1
> 3  -32 -22   1
> 4  -33 -22   1
> 5  -32 -23   2
> 6  -31 -23   2
> 7  -31 -22   2
> 8  -32 -22   2
> 9  -31 -23   3
> 10 -30 -23   3
> 11 -30 -22   3
> 12 -31 -22   3
> 13 -33 -22   4
> 14 -32 -22   4
> 15 -32 -21   4
> 16 -33 -21   4
> 17 -32 -22   5
> 18 -31 -22   5
> 19 -31 -21   5
> 20 -32 -21   5
> 21 -31 -22   6
> 22 -30 -22   6
> 23 -30 -21   6
> 24 -31 -21   6
>
>         [[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

how to create several polygons from a list of vertices

Tue, 08/14/2018 - 15:16
Hi,

I have a data.frame with the vertices (lon / lat) and codes from several
squares (more than 500 in the real dataset).
I want to create an object with these polygons (squares) and after this
export it as a shapefile.
With the script below I can draw one square.
library(sp)
P1 = Polygon(vertices[1:4,1:2])
Ps1 = SpatialPolygons(list(Polygons(list(P1), ID = "1")),
proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
plot(Ps1, axes = TRUE)

Now I'm trying to create one object with all squares at once.
Is it possible?

Thanks a lot,

Antônio Olinto

sample data:vertices
   lon lat cod
1  -33 -23   1
2  -32 -23   1
3  -32 -22   1
4  -33 -22   1
5  -32 -23   2
6  -31 -23   2
7  -31 -22   2
8  -32 -22   2
9  -31 -23   3
10 -30 -23   3
11 -30 -22   3
12 -31 -22   3
13 -33 -22   4
14 -32 -22   4
15 -32 -21   4
16 -33 -21   4
17 -32 -22   5
18 -31 -22   5
19 -31 -21   5
20 -32 -21   5
21 -31 -22   6
22 -30 -22   6
23 -30 -21   6
24 -31 -21   6

        [[alternative HTML version deleted]]

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

Re: Extract coordinates from rasterbrick

Tue, 08/14/2018 - 14:32
Dear Ákos and Vijay,

Thank you very much for the solutions. They work perfectly!

Sincerely,

Milu

On Tue, Aug 14, 2018 at 6:12 PM, Vijay Lulla <[hidden email]> wrote:

> And if you need coordinates as part of your data frame just do
> cbind(coordinates(x), as.data.frame(x, xy = TRUE))
>
> On Tue, Aug 14, 2018 at 11:50 AM Bede-Fazekas Ákos <[hidden email]>
> wrote:
>
> > Dear Milu,
> >
> > I think that you are looking for as.data.frame(x, xy = TRUE).
> >
> > HTH,
> > Ákos Bede-Fazekas
> > Hungarian Academy of Sciences
> >
> >
> > 2018.08.14. 17:03 keltezéssel, Miluji Sb írta:
> > > Dear all,
> > >
> > > I have the following rasterbrick (x);
> > >
> > > class       : RasterBrick
> > > dimensions  : 112, 272, 30464, 7305  (nrow, ncol, ncell, nlayers)
> > > resolution  : 0.25, 0.25  (x, y)
> > > extent      : -132, -64, 24, 52  (xmin, xmax, ymin, ymax)
> > > coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> > > data source : /projectnb/climpct/TEX_1986_2005.nc
> > > names       : X1986.01.01, X1986.01.02, X1986.01.03, X1986.01.04,
> > > X1986.01.05, X1986.01.06, X1986.01.07, X1986.01.08, X1986.01.09,
> > > X1986.01.10, X1986.01.11, X1986.01.12, X1986.01.13, X1986.01.14,
> > > X1986.01.15, ...
> > > Date        : 1986-01-01, 2005-12-31 (min, max)
> > > varname     : tex
> > >
> > >
> > > I can obtain the values by;
> > >
> > > as.data.frame(getValues(x))
> > >
> > > However, I would like to extract the values by all coordinates in the
> > file
> > > in the form of a dataframe. Is it possible to do so? Any help will be
> > > greatly appreciated. Thank you.
> > >
> > > Sincerely,
> > >
> > > Milu
> > >
> > >       [[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
> >
>
>         [[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: Extract coordinates from rasterbrick

Tue, 08/14/2018 - 11:12
And if you need coordinates as part of your data frame just do
cbind(coordinates(x), as.data.frame(x, xy = TRUE))

On Tue, Aug 14, 2018 at 11:50 AM Bede-Fazekas Ákos <[hidden email]>
wrote:

> Dear Milu,
>
> I think that you are looking for as.data.frame(x, xy = TRUE).
>
> HTH,
> Ákos Bede-Fazekas
> Hungarian Academy of Sciences
>
>
> 2018.08.14. 17:03 keltezéssel, Miluji Sb írta:
> > Dear all,
> >
> > I have the following rasterbrick (x);
> >
> > class       : RasterBrick
> > dimensions  : 112, 272, 30464, 7305  (nrow, ncol, ncell, nlayers)
> > resolution  : 0.25, 0.25  (x, y)
> > extent      : -132, -64, 24, 52  (xmin, xmax, ymin, ymax)
> > coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> > data source : /projectnb/climpct/TEX_1986_2005.nc
> > names       : X1986.01.01, X1986.01.02, X1986.01.03, X1986.01.04,
> > X1986.01.05, X1986.01.06, X1986.01.07, X1986.01.08, X1986.01.09,
> > X1986.01.10, X1986.01.11, X1986.01.12, X1986.01.13, X1986.01.14,
> > X1986.01.15, ...
> > Date        : 1986-01-01, 2005-12-31 (min, max)
> > varname     : tex
> >
> >
> > I can obtain the values by;
> >
> > as.data.frame(getValues(x))
> >
> > However, I would like to extract the values by all coordinates in the
> file
> > in the form of a dataframe. Is it possible to do so? Any help will be
> > greatly appreciated. Thank you.
> >
> > Sincerely,
> >
> > Milu
> >
> >       [[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
>
        [[alternative HTML version deleted]]

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

Re: Extract coordinates from rasterbrick

Tue, 08/14/2018 - 10:49
Dear Milu,

I think that you are looking for as.data.frame(x, xy = TRUE).

HTH,
Ákos Bede-Fazekas
Hungarian Academy of Sciences


2018.08.14. 17:03 keltezéssel, Miluji Sb írta:
> Dear all,
>
> I have the following rasterbrick (x);
>
> class       : RasterBrick
> dimensions  : 112, 272, 30464, 7305  (nrow, ncol, ncell, nlayers)
> resolution  : 0.25, 0.25  (x, y)
> extent      : -132, -64, 24, 52  (xmin, xmax, ymin, ymax)
> coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> data source : /projectnb/climpct/TEX_1986_2005.nc
> names       : X1986.01.01, X1986.01.02, X1986.01.03, X1986.01.04,
> X1986.01.05, X1986.01.06, X1986.01.07, X1986.01.08, X1986.01.09,
> X1986.01.10, X1986.01.11, X1986.01.12, X1986.01.13, X1986.01.14,
> X1986.01.15, ...
> Date        : 1986-01-01, 2005-12-31 (min, max)
> varname     : tex
>
>
> I can obtain the values by;
>
> as.data.frame(getValues(x))
>
> However, I would like to extract the values by all coordinates in the file
> in the form of a dataframe. Is it possible to do so? Any help will be
> greatly appreciated. Thank you.
>
> Sincerely,
>
> Milu
>
> [[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

Extract coordinates from rasterbrick

Tue, 08/14/2018 - 10:03
Dear all,

I have the following rasterbrick (x);

class       : RasterBrick
dimensions  : 112, 272, 30464, 7305  (nrow, ncol, ncell, nlayers)
resolution  : 0.25, 0.25  (x, y)
extent      : -132, -64, 24, 52  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : /projectnb/climpct/TEX_1986_2005.nc
names       : X1986.01.01, X1986.01.02, X1986.01.03, X1986.01.04,
X1986.01.05, X1986.01.06, X1986.01.07, X1986.01.08, X1986.01.09,
X1986.01.10, X1986.01.11, X1986.01.12, X1986.01.13, X1986.01.14,
X1986.01.15, ...
Date        : 1986-01-01, 2005-12-31 (min, max)
varname     : tex


I can obtain the values by;

as.data.frame(getValues(x))

However, I would like to extract the values by all coordinates in the file
in the form of a dataframe. Is it possible to do so? Any help will be
greatly appreciated. Thank you.

Sincerely,

Milu

        [[alternative HTML version deleted]]

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

Re: How do i plot a graph with x axis vertically labelled?

Mon, 08/13/2018 - 09:32
The graphics parameters are all in the help file for par
?par

The one you need is las.

Without trying to wade thru your mangled HTML data (dput() is the best
way to provide reproducible data), check this out:

barplot(runif(5), names.arg = c("25/08/2004", "26/09/2004",
"02/11/2004", "11/11/2004", "22/11/2004"), main = 'Females 2005', ylab
= ' Age/Days', xlab = 'Birth/Hatch date', las=2)

You'll need to mess with label positioning, and with margin size, to
get everything looking as you wish.

But I'm not sure this is the most effective way to show
irregularly-spaced dates. Why not convert your dates to Date or other
time series, and plot them to scale as points?

Sarah
On Mon, Aug 13, 2018 at 6:25 AM Vasana Tutjavi via R-sig-Geo
<[hidden email]> wrote:
>
> Dear R-sig-geo experts,
>
> I want to plot a bar graph with the labels on the x axis in a vertical format. The data is the Ages of Squid (on the Y axis) and the date of Birth on the X axis, I am struggling to change the labels on the x axis to a vertical format and I want all the dates to show at all the bars. My data is as follows:
> Birth
> Sex
> Age
> 25/08/2004
> F
> 2
> 26/09/2004
> F
> 3
> 02/11/2004
> F
> 3
> 11/11/2004
> F
> 1
> 22/11/2004
> F
> 3
> 04/12/2004
> F
> 1
> 08/12/2004
> F
> 4
> 15/12/2004
> F
> 4
> 21/12/2004
> F
> 5
> 24/12/2004
> F
> 1
> 04/01/2005
> F
> 1
> 06/01/2005
> F
> 3
> 13/01/2005
> F
> 2
> 18/01/2005
> F
> 1
> 19/01/2005
> F
> 3
> 09/02/2005
> F
> 2
> 11/02/2005
> F
> 1
> 13/02/2005
> F
> 1
> 17/02/2005
> F
> 1
> 18/02/2005
> F
> 8
> 19/02/2005
> F
> 5
> 24/02/2005
> F
> 1
> 28/02/2005
> F
> 42
> 01/03/2005
> F
> 1
> 07/03/2005
> F
> 1
> 08/03/2005
> F
> 2
> 09/03/2005
> F
> 1
> 14/03/2005
> F
> 7
> 16/03/2005
> F
> 1
> 25/03/2005
> F
> 3
> 01/04/2005
> F
> 1
> 06/04/2005
> F
> 5
> 07/04/2005
> F
> 1
> 08/04/2005
> F
> 1
> 13/04/2005
> F
> 1
> 14/04/2005
> F
> 1
> 27/04/2005
> F
> 1
> 01/05/2005
> F
> 2
> 07/05/2005
> F
> 1
> 08/05/2005
> F
> 2
> 16/05/2005
> F
> 1
> 30/05/2005
> F
> 8
> 06/06/2005
> F
> 1
> Here is the script I tried:
>
> barplot(Females2005$Age,  names.arg = c("25/08/2004", "26/09/2004", "02/11/2004", "11/11/2004", "22/11/2004", "04/12/2004", "08/12/2004",
>                                        "15/12/2004", "21/12/2004", "24/12/2004", "04/01/2005", "06/01/2005", "13/01/2005", "18/01/2005",
>                                        "19/01/2005", "09/02/2005", "11/02/2005", "13/02/2005", "17/02/2005", "18/02/2005", "19/02/2005",
>                                       "24/02/2005", "28/02/2005", "01/03/2005", "07/03/2005", "08/03/2005", "09/03/2005", "14/03/2005",
>                                     "16/03/2005", "25/03/2005", "01/04/2005", "06/04/2005", "07/04/2005", "08/04/2005", "13/04/2005",
>                                       "14/04/2005", "27/04/2005", "01/05/2005", "07/05/2005", "08/05/2005", "16/05/2005", "30/05/2005",
>                                     "06/06/2005"), main = 'Females 2005', ylab = 'Age/Days', xlab = 'Birth/Hatch date')
>
> Many thanks
> Vasana
> Sent from Mail for Windows 10
>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo


--
Sarah Goslee
http://www.functionaldiversity.org

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

Pages