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 38 min ago

Re: Geographically weighted regression

Fri, 02/22/2019 - 12:50
The first step should be to look at

str(Daten90)
str(Daten10)

and if that doesn't solve the problem, then consider a reproducible
example, or at the very least posting the results of the above to this
list.

Sarah

On Fri, Feb 22, 2019 at 7:38 AM <[hidden email]> wrote:
>
> Dear all,
>
> I am currently working out a geographically weighted regression, in which 90% of the data set the model should be calculated and for 10% of the values to be predicted. For the prediction I use the function gwr.predict from the package GWModel:
>
>  Erg<-gwr.predict(formula=Ziel~ as.factor(Var1) + log(Var2, base = exp(1)) + Var3, data = Daten90,predictdata = Daten10,bw = bwG, kernel = "gaussian",adaptive = FALSE, p = 2, theta = 0, longlat = FALSE)
>
> I always get this error, although Daten10 and Daten90 have the same structure:
> Error in gwr.predict(formula = Ziel~ as.factor(Var1) + log(Var2, base = exp(1)) + Var3, :
> All the independent variables should be included in the predictdata.
>
> Can you tell me what the problem with this code is?
> Or is there any other way for a GWR and the prediction?
>
> Thank you,
> Christoph
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo


--
Sarah Goslee (she/her)
http://www.numberwright.com

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

Re: create a raster time series with negative years (years Before Present)

Fri, 02/22/2019 - 09:12
See ?as.Date:

"Years before 1CE (aka 1AD) will probably not be handled correctly."

On 2/22/19 3:59 PM, Jackson Rodrigues wrote:
> Hey guys,
>
> I need some help to set a time serie (ts) to raster stack, however my
> raster serie has negative ages (years Before Present (yrBP)).
>
> I have tried several codes and packages unsuccessfully.
>
> If succeed, I would like to apply seasonal and other calculations to my
> new  raster time serie.
>
> Below are some codes from package rts. I adapted years to negative value
> only as a example.
>
> #####
>> path <- system.file("external", package="rts") # location of files
>> lst <- list.files(path=path,pattern='.asc$',full.names=TRUE)
>> lst # list of raster files
>> r <- stack(lst) # creating a RasterStack object
>> r
>> d <- c("-2000-02-01","-2000-03-01","-2000-04-01","-2000-05-01") #
> corresponding dates to 4 rasters
>> d <- as.Date(d) # or d <- as.POSIXct(d)
> Error in charToDate(x) :
>   character string is not in a standard unambiguous format
>> rt <- rts(r,d) # creating a RasterStackTS object
>> rt
>> plot(rt)
> #########
>
> best,
>
> Jackson Rodrigues
> --
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> --
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081

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

pEpkey.asc (2K) Download Attachment

create a raster time series with negative years (years Before Present)

Fri, 02/22/2019 - 08:59
Hey guys,

I need some help to set a time serie (ts) to raster stack, however my
raster serie has negative ages (years Before Present (yrBP)).

I have tried several codes and packages unsuccessfully.

If succeed, I would like to apply seasonal and other calculations to my
new  raster time serie.

Below are some codes from package rts. I adapted years to negative value
only as a example.

#####
>path <- system.file("external", package="rts") # location of files
>lst <- list.files(path=path,pattern='.asc$',full.names=TRUE)
>lst # list of raster files
>r <- stack(lst) # creating a RasterStack object
>r
>d <- c("-2000-02-01","-2000-03-01","-2000-04-01","-2000-05-01") #
corresponding dates to 4 rasters
>d <- as.Date(d) # or d <- as.POSIXct(d)
Error in charToDate(x) :
  character string is not in a standard unambiguous format
>rt <- rts(r,d) # creating a RasterStackTS object
>rt
>plot(rt)
#########

best,

Jackson Rodrigues
--

        [[alternative HTML version deleted]]

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

Geographically weighted regression

Fri, 02/22/2019 - 06:37
Dear all,
  I am currently working out a geographically weighted regression, in which 90% of the data set the model should be calculated and for 10% of the values to be predicted. For the prediction I use the function gwr.predict from the package GWModel:
 Erg<-gwr.predict(formula=Ziel~ as.factor(Var1) + log(Var2, base = exp(1)) + Var3, data = Daten90,predictdata = Daten10,bw = bwG, kernel = "gaussian",adaptive = FALSE, p = 2, theta = 0, longlat = FALSE)

I always get this error, although Daten10 and Daten90 have the same structure:
Error in gwr.predict(formula = Ziel~ as.factor(Var1) + log(Var2, base = exp(1)) + Var3, : All the independent variables should be included in the predictdata.   Can you tell me what the problem with this code is?
Or is there any other way for a GWR and the prediction?

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

Re: Use rbind for more than two data frame

Thu, 02/21/2019 - 06:55
Thank you!
Sorry for my silly doubt.

Peter van Horssen <[hidden email]> escreveu no dia quinta,
21/02/2019 à(s) 11:39:

> Yes,
>
> alter columname   "bio15_FP" to "bio15" in  table 'sdmdata_F'
>
> :-)
>
> greetings  Peter
>
>
> Op 21-2-2019 om 13:29 schreef Lara Silva:
> > Hello,
> >
> > I am trying to join different data frames using rbind but I have the
> > following mensage:
> >
> > sdmata_global <-
> > rbind(sdmdata_SM,sdmdata_P,sdmdata_T,sdmdata_F,sdmdata_SJ)
> >
> > Error in match.names(clabs, names(xi)) :
> > names do not match previous names
> >
> >> class(sdmdata_SM)
> > [1] "data.frame"
> >> class(sdmdata_P)
> > [1] "data.frame"
> >> class(sdmdata_T)
> > [1] "data.frame"
> >> class(sdmdata_F)
> > [1] "data.frame"
> >> class(sdmdata_SJ)
> > [1] "data.frame"
> >> dim(sdmdata_SM)
> > [1] 10003    22
> >> dim(sdmdata_P)
> > [1] 10005    22
> >> dim(sdmdata_T)
> > [1] 10004    22
> >> dim(sdmdata_F)
> > [1] 10001    22
> >> dim(sdmdata_SJ)
> > [1] 10001    22
> >
> >
> >> dim(sdmdata_SM)
> > [1] 10003    22
> >> dim(sdmdata_SM)
> > [1] 10003    22
> >> dim(sdmdata_P)
> > [1] 10005    22
> >> dim(sdmdata_T)
> > [1] 10004    22
> >> dim(sdmdata_F)
> > [1] 10001    22
> >> dim(sdmdata_SJ)
> > [1] 10001    22
> >> names(sdmdata_SM)
> >   [1] "pb"    "bio1"  "bio10" "bio11" "bio12" "bio13" "bio14" "bio15"
> "bio16"
> > [10] "bio17" "bio18" "bio19" "bio2"  "bio3"  "bio4"  "bio4a" "bio5"
> "bio6"
> > [19] "bio7"  "bio8"  "bio9"  "code"
> >> names(sdmdata_P)
> >   [1] "pb"    "bio1"  "bio10" "bio11" "bio12" "bio13" "bio14" "bio15"
> "bio16"
> > [10] "bio17" "bio18" "bio19" "bio2"  "bio3"  "bio4"  "bio4a" "bio5"
> "bio6"
> > [19] "bio7"  "bio8"  "bio9"  "code"
> >> names(sdmdata_T)
> >   [1] "pb"    "bio1"  "bio10" "bio11" "bio12" "bio13" "bio14" "bio15"
> "bio16"
> > [10] "bio17" "bio18" "bio19" "bio2"  "bio3"  "bio4"  "bio4a" "bio5"
> "bio6"
> > [19] "bio7"  "bio8"  "bio9"  "code"
> >> names(sdmdata_F)
> >   [1] "pb"       "bio1"     "bio10"    "bio11"    "bio12"    "bio13"
> >   [7] "bio14"    "bio15_FP" "bio16"    "bio17"    "bio18"    "bio19"
> > [13] "bio2"     "bio3"     "bio4"     "bio4a"    "bio5"     "bio6"
> > [19] "bio7"     "bio8"     "bio9"     "code"
> >> names(sdmdata_SJ)
> >   [1] "pb"    "bio1"  "bio10" "bio11" "bio12" "bio13" "bio14" "bio15"
> "bio16"
> > [10] "bio17" "bio18" "bio19" "bio2"  "bio3"  "bio4"  "bio4a" "bio5"
> "bio6"
> > [19] "bio7"  "bio8"  "bio9"  "code"
> > Any suggestion?
> >
> > Thanks!
> >
> > Lara
> >
> >       [[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: Use rbind for more than two data frame

Thu, 02/21/2019 - 06:37
Yes,

alter columname   "bio15_FP" to "bio15" in  table 'sdmdata_F'

:-)

greetings  Peter


Op 21-2-2019 om 13:29 schreef Lara Silva:
> Hello,
>
> I am trying to join different data frames using rbind but I have the
> following mensage:
>
> sdmata_global <-
> rbind(sdmdata_SM,sdmdata_P,sdmdata_T,sdmdata_F,sdmdata_SJ)
>
> Error in match.names(clabs, names(xi)) :
> names do not match previous names
>
>> class(sdmdata_SM)
> [1] "data.frame"
>> class(sdmdata_P)
> [1] "data.frame"
>> class(sdmdata_T)
> [1] "data.frame"
>> class(sdmdata_F)
> [1] "data.frame"
>> class(sdmdata_SJ)
> [1] "data.frame"
>> dim(sdmdata_SM)
> [1] 10003    22
>> dim(sdmdata_P)
> [1] 10005    22
>> dim(sdmdata_T)
> [1] 10004    22
>> dim(sdmdata_F)
> [1] 10001    22
>> dim(sdmdata_SJ)
> [1] 10001    22
>
>
>> dim(sdmdata_SM)
> [1] 10003    22
>> dim(sdmdata_SM)
> [1] 10003    22
>> dim(sdmdata_P)
> [1] 10005    22
>> dim(sdmdata_T)
> [1] 10004    22
>> dim(sdmdata_F)
> [1] 10001    22
>> dim(sdmdata_SJ)
> [1] 10001    22
>> names(sdmdata_SM)
>   [1] "pb"    "bio1"  "bio10" "bio11" "bio12" "bio13" "bio14" "bio15" "bio16"
> [10] "bio17" "bio18" "bio19" "bio2"  "bio3"  "bio4"  "bio4a" "bio5"  "bio6"
> [19] "bio7"  "bio8"  "bio9"  "code"
>> names(sdmdata_P)
>   [1] "pb"    "bio1"  "bio10" "bio11" "bio12" "bio13" "bio14" "bio15" "bio16"
> [10] "bio17" "bio18" "bio19" "bio2"  "bio3"  "bio4"  "bio4a" "bio5"  "bio6"
> [19] "bio7"  "bio8"  "bio9"  "code"
>> names(sdmdata_T)
>   [1] "pb"    "bio1"  "bio10" "bio11" "bio12" "bio13" "bio14" "bio15" "bio16"
> [10] "bio17" "bio18" "bio19" "bio2"  "bio3"  "bio4"  "bio4a" "bio5"  "bio6"
> [19] "bio7"  "bio8"  "bio9"  "code"
>> names(sdmdata_F)
>   [1] "pb"       "bio1"     "bio10"    "bio11"    "bio12"    "bio13"
>   [7] "bio14"    "bio15_FP" "bio16"    "bio17"    "bio18"    "bio19"
> [13] "bio2"     "bio3"     "bio4"     "bio4a"    "bio5"     "bio6"
> [19] "bio7"     "bio8"     "bio9"     "code"
>> names(sdmdata_SJ)
>   [1] "pb"    "bio1"  "bio10" "bio11" "bio12" "bio13" "bio14" "bio15" "bio16"
> [10] "bio17" "bio18" "bio19" "bio2"  "bio3"  "bio4"  "bio4a" "bio5"  "bio6"
> [19] "bio7"  "bio8"  "bio9"  "code"
> Any suggestion?
>
> Thanks!
>
> Lara
>
> [[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

Use rbind for more than two data frame

Thu, 02/21/2019 - 06:29
Hello,

I am trying to join different data frames using rbind but I have the
following mensage:

sdmata_global <-
rbind(sdmdata_SM,sdmdata_P,sdmdata_T,sdmdata_F,sdmdata_SJ)

Error in match.names(clabs, names(xi)) :
names do not match previous names

> class(sdmdata_SM)
[1] "data.frame"
> class(sdmdata_P)
[1] "data.frame"
> class(sdmdata_T)
[1] "data.frame"
> class(sdmdata_F)
[1] "data.frame"
> class(sdmdata_SJ)
[1] "data.frame"
> dim(sdmdata_SM)
[1] 10003    22
> dim(sdmdata_P)
[1] 10005    22
> dim(sdmdata_T)
[1] 10004    22
> dim(sdmdata_F)
[1] 10001    22
> dim(sdmdata_SJ)
[1] 10001    22


> dim(sdmdata_SM)
[1] 10003    22
> dim(sdmdata_SM)
[1] 10003    22
> dim(sdmdata_P)
[1] 10005    22
> dim(sdmdata_T)
[1] 10004    22
> dim(sdmdata_F)
[1] 10001    22
> dim(sdmdata_SJ)
[1] 10001    22
>

> names(sdmdata_SM)
 [1] "pb"    "bio1"  "bio10" "bio11" "bio12" "bio13" "bio14" "bio15" "bio16"
[10] "bio17" "bio18" "bio19" "bio2"  "bio3"  "bio4"  "bio4a" "bio5"  "bio6"
[19] "bio7"  "bio8"  "bio9"  "code"
> names(sdmdata_P)
 [1] "pb"    "bio1"  "bio10" "bio11" "bio12" "bio13" "bio14" "bio15" "bio16"
[10] "bio17" "bio18" "bio19" "bio2"  "bio3"  "bio4"  "bio4a" "bio5"  "bio6"
[19] "bio7"  "bio8"  "bio9"  "code"
> names(sdmdata_T)
 [1] "pb"    "bio1"  "bio10" "bio11" "bio12" "bio13" "bio14" "bio15" "bio16"
[10] "bio17" "bio18" "bio19" "bio2"  "bio3"  "bio4"  "bio4a" "bio5"  "bio6"
[19] "bio7"  "bio8"  "bio9"  "code"
> names(sdmdata_F)
 [1] "pb"       "bio1"     "bio10"    "bio11"    "bio12"    "bio13"
 [7] "bio14"    "bio15_FP" "bio16"    "bio17"    "bio18"    "bio19"
[13] "bio2"     "bio3"     "bio4"     "bio4a"    "bio5"     "bio6"
[19] "bio7"     "bio8"     "bio9"     "code"
> names(sdmdata_SJ)
 [1] "pb"    "bio1"  "bio10" "bio11" "bio12" "bio13" "bio14" "bio15" "bio16"
[10] "bio17" "bio18" "bio19" "bio2"  "bio3"  "bio4"  "bio4a" "bio5"  "bio6"
[19] "bio7"  "bio8"  "bio9"  "code"
>

Any suggestion?

Thanks!

Lara

        [[alternative HTML version deleted]]

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

Re: Reverse color in stplot

Thu, 02/21/2019 - 01:45
I think that it accepts an argument col.regions where you specify a
color palette; if you have a particular one, called pal, and want to
revert it, use rev(pal).

On 2/21/19 7:02 AM, Mahsa Sameti wrote:
> Dear all
>
> I want to reverse color in my spatio-temporal maps that are created with
> stplot function. How can i solve this problem?
>
> Thanks alot
> Mahsa Sameti
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> --
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081

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

pEpkey.asc (2K) Download Attachment

Reverse color in stplot

Wed, 02/20/2019 - 23:54
Dear all

I want to reverse color in my spatio-temporal maps that are created with
stplot function. How can i solve this problem?

Thanks alot
Mahsa Sameti

        [[alternative HTML version deleted]]

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

Re: [R-sig-eco] Problem in generate random samples in r

Wed, 02/20/2019 - 05:29
Thanks!

Lara

Henrik Eckermann <[hidden email]> escreveu no dia quarta,
20/02/2019 à(s) 10:25:

> Hi Lara,
>
> the error is quite informative. You use the sample function twice. At
> least for one the length of the vector where you sample from is shorter
> (contains less values) than than the number of samples you wanna draw.
>
> best,
>
> Henrik
>
> > On 20. Feb 2019, at 12:20, Lara Silva <[hidden email]> wrote:
> >
> > Hello,
> >
> > I am trying to generate random samples from the following code
> >
> > ### Setting random seed to always create the same
> > ### Random set of points
> > set.seed(0)
> >
> > absences_15000<-absences[sample(nrow(absences), 15000),]
> > points(absences_15000, cex=0.1)
> >
> > ## Subsample_10000
> > set.seed(0)
> >
> > absences_10000<-absences_15000[sample(nrow(absences_15000), 10000),]
> > dim(absences_10000)
> >
> > I get the following error:
> >
> > Error in sample.int(length(x), size, replace, prob) :
> >  cannot take a sample larger than the population when 'replace = FALSE'
> >
> > Any advice?
> >
> > Regards,
> >
> > Lara
> >
> >       [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-ecology mailing list
> > [hidden email]
> > https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>
>
        [[alternative HTML version deleted]]

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

Re: [R-sig-eco] Problem in generate random samples in r

Wed, 02/20/2019 - 05:28
My guess would be that absences_150000 doesn't have 10000 rows. Can you
confirm or refute this?

Cheers,
Roman

On Wed, Feb 20, 2019 at 12:20 PM Lara Silva <[hidden email]>
wrote:

> Hello,
>
> I am trying to generate random samples from the following code
>
> ### Setting random seed to always create the same
> ### Random set of points
> set.seed(0)
>
> absences_15000<-absences[sample(nrow(absences), 15000),]
> points(absences_15000, cex=0.1)
>
> ## Subsample_10000
> set.seed(0)
>
> absences_10000<-absences_15000[sample(nrow(absences_15000), 10000),]
> dim(absences_10000)
>
> I get the following error:
>
> Error in sample.int(length(x), size, replace, prob) :
>   cannot take a sample larger than the population when 'replace = FALSE'
>
> Any advice?
>
> Regards,
>
> Lara
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-ecology mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>

--
In God we trust, all others bring data.

        [[alternative HTML version deleted]]

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

Re: Plot coordinates on world map with Robinson CRS - ggplot2

Mon, 02/18/2019 - 05:43
Thank you very much. This works perfectly. Appreciate it.

Sincerely,

Milu

On Mon, Feb 18, 2019 at 11:01 AM Michael Sumner <[hidden email]> wrote:

> You are very close, just needed the expert-knowledge for the incantation
> required to transform the longlat points (as was already done for the
> prj.coord, btw).
>
> HTH: https://gist.github.com/mdsumner/afd9228aa1dc4652f8bf193cad4245aa
>
> It does come as a surprise to many that R's plotting is not unit or
> projection aware (well, it is sometimes, but not usually ... a modern
> example is ggplot2::geom_sf  and the sf package but you need quite a lot of
> prior-won expertise there too, good luck).
>
> Cheers, Mike.
>
> On Mon, 18 Feb 2019 at 20:01 Miluji Sb <[hidden email]> wrote:
>
>> Dear all,
>>
>> I am trying to plot coordinates on a world map with Robinson CRS. While
>> the
>> world map is generated without any issues, when I try to plot the points -
>> I only get a single point.
>>
>> The code I am using and the coordinates data is below. What am I doing
>> wrong? Any help/suggestions will be highly appreciated.
>>
>> library(data.table)
>> library(foreign)
>> library(readstata13)
>> library(rgdal)
>> library(maptools)
>> library(ggplot2)
>> library(dplyr)
>>
>> load(url("
>>
>> https://github.com/valentinitnelav/RandomScripts/blob/master/NaturalEarth.RData?raw=true
>> "))
>>
>> PROJ <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84
>> +units=m +no_defs"
>>
>> NE_countries_rob  <- spTransform(NE_countries, CRSobj = PROJ)
>> NE_graticules_rob <- spTransform(NE_graticules, CRSobj = PROJ)
>> NE_box_rob        <- spTransform(NE_box, CRSobj = PROJ)
>>
>> # project long-lat coordinates for graticule label data frames
>> # (two extra columns with projected XY are created)
>> prj.coord <- project(cbind(lbl.Y$lon, lbl.Y$lat), proj=PROJ)
>> lbl.Y.prj <- cbind(prj.coord, lbl.Y)
>> names(lbl.Y.prj)[1:2] <- c("X.prj","Y.prj")
>>
>> prj.coord <- project(cbind(lbl.X$lon, lbl.X$lat), proj=PROJ)
>> lbl.X.prj <- cbind(prj.coord, lbl.X)
>> names(lbl.X.prj)[1:2] <- c("X.prj","Y.prj")
>>
>> m <- ggplot() +
>>   # add Natural Earth countries projected to Robinson, give black border
>> and fill with gray
>>   geom_polygon(data=NE_countries_rob, aes(long,lat, group=group),
>> colour="black", fill="gray80", size = 0.25) +
>>   # Note: "Regions defined for each Polygons" warning has to do with
>> fortify transformation. Might get deprecated in future!
>>   # alternatively, use use map_data(NE_countries) to transform to data
>> frame and then use project() to change to desired projection.
>>   # add Natural Earth box projected to Robinson
>>   geom_polygon(data=NE_box_rob, aes(x=long, y=lat), colour="black",
>> fill="transparent", size = 0.25) +
>>   # add graticules projected to Robinson
>>   geom_path(data=NE_graticules_rob, aes(long, lat, group=group),
>> linetype="dotted", color="grey50", size = 0.25) +
>>   # add graticule labels - latitude and longitude
>>   geom_text(data = lbl.Y.prj, aes(x = X.prj, y = Y.prj, label = lbl),
>> color="grey50", size=2) +
>>   geom_text(data = lbl.X.prj, aes(x = X.prj, y = Y.prj, label = lbl),
>> color="grey50", size=2) +
>>   # the default, ratio = 1 in coord_fixed ensures that one unit on the
>> x-axis is the same length as one unit on the y-axis
>>   coord_fixed(ratio = 1) +
>>   geom_point(data=df,
>>              aes(x=lon, y=lat), colour="Deep Pink",
>>              fill="Pink",pch=21, size=2, alpha=I(0.7))
>>   # remove the background and default gridlines
>>   theme_void()
>>
>> ## coordinates dataframe
>> dput(df)
>> structure(list(lon = c(2.67569724621467, 17.5766416259819,
>> 28.4126232192772,
>> 23.8147674538232, 29.8917589327105), lat = c(28.1503115976162,
>> -12.3388787380201, 9.78891068739477, -22.1873831176644, -3.36546931479253
>> )), class = "data.frame", row.names = c(NA, -5L))
>>
>> Sincerely,
>>
>> Milu
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
> --
> Dr. Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> 203 Channel Highway
> Kingston Tasmania 7050 Australia
>
>
        [[alternative HTML version deleted]]

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

Re: Plot coordinates on world map with Robinson CRS - ggplot2

Mon, 02/18/2019 - 04:00
You are very close, just needed the expert-knowledge for the incantation
required to transform the longlat points (as was already done for the
prj.coord, btw).

HTH: https://gist.github.com/mdsumner/afd9228aa1dc4652f8bf193cad4245aa

It does come as a surprise to many that R's plotting is not unit or
projection aware (well, it is sometimes, but not usually ... a modern
example is ggplot2::geom_sf  and the sf package but you need quite a lot of
prior-won expertise there too, good luck).

Cheers, Mike.

On Mon, 18 Feb 2019 at 20:01 Miluji Sb <[hidden email]> wrote:

> Dear all,
>
> I am trying to plot coordinates on a world map with Robinson CRS. While the
> world map is generated without any issues, when I try to plot the points -
> I only get a single point.
>
> The code I am using and the coordinates data is below. What am I doing
> wrong? Any help/suggestions will be highly appreciated.
>
> library(data.table)
> library(foreign)
> library(readstata13)
> library(rgdal)
> library(maptools)
> library(ggplot2)
> library(dplyr)
>
> load(url("
>
> https://github.com/valentinitnelav/RandomScripts/blob/master/NaturalEarth.RData?raw=true
> "))
>
> PROJ <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84
> +units=m +no_defs"
>
> NE_countries_rob  <- spTransform(NE_countries, CRSobj = PROJ)
> NE_graticules_rob <- spTransform(NE_graticules, CRSobj = PROJ)
> NE_box_rob        <- spTransform(NE_box, CRSobj = PROJ)
>
> # project long-lat coordinates for graticule label data frames
> # (two extra columns with projected XY are created)
> prj.coord <- project(cbind(lbl.Y$lon, lbl.Y$lat), proj=PROJ)
> lbl.Y.prj <- cbind(prj.coord, lbl.Y)
> names(lbl.Y.prj)[1:2] <- c("X.prj","Y.prj")
>
> prj.coord <- project(cbind(lbl.X$lon, lbl.X$lat), proj=PROJ)
> lbl.X.prj <- cbind(prj.coord, lbl.X)
> names(lbl.X.prj)[1:2] <- c("X.prj","Y.prj")
>
> m <- ggplot() +
>   # add Natural Earth countries projected to Robinson, give black border
> and fill with gray
>   geom_polygon(data=NE_countries_rob, aes(long,lat, group=group),
> colour="black", fill="gray80", size = 0.25) +
>   # Note: "Regions defined for each Polygons" warning has to do with
> fortify transformation. Might get deprecated in future!
>   # alternatively, use use map_data(NE_countries) to transform to data
> frame and then use project() to change to desired projection.
>   # add Natural Earth box projected to Robinson
>   geom_polygon(data=NE_box_rob, aes(x=long, y=lat), colour="black",
> fill="transparent", size = 0.25) +
>   # add graticules projected to Robinson
>   geom_path(data=NE_graticules_rob, aes(long, lat, group=group),
> linetype="dotted", color="grey50", size = 0.25) +
>   # add graticule labels - latitude and longitude
>   geom_text(data = lbl.Y.prj, aes(x = X.prj, y = Y.prj, label = lbl),
> color="grey50", size=2) +
>   geom_text(data = lbl.X.prj, aes(x = X.prj, y = Y.prj, label = lbl),
> color="grey50", size=2) +
>   # the default, ratio = 1 in coord_fixed ensures that one unit on the
> x-axis is the same length as one unit on the y-axis
>   coord_fixed(ratio = 1) +
>   geom_point(data=df,
>              aes(x=lon, y=lat), colour="Deep Pink",
>              fill="Pink",pch=21, size=2, alpha=I(0.7))
>   # remove the background and default gridlines
>   theme_void()
>
> ## coordinates dataframe
> dput(df)
> structure(list(lon = c(2.67569724621467, 17.5766416259819,
> 28.4126232192772,
> 23.8147674538232, 29.8917589327105), lat = c(28.1503115976162,
> -12.3388787380201, 9.78891068739477, -22.1873831176644, -3.36546931479253
> )), class = "data.frame", row.names = c(NA, -5L))
>
> Sincerely,
>
> Milu
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> --
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

        [[alternative HTML version deleted]]

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

Re: Plot coordinates on world map with Robinson CRS - ggplot2

Mon, 02/18/2019 - 03:17
On Mon, 18 Feb 2019, Miluji Sb wrote:

> Dear all,
>
> I am trying to plot coordinates on a world map with Robinson CRS. While the
> world map is generated without any issues, when I try to plot the points -
> I only get a single point.

Well, the point is correct, since you didn't project df, did you? This
seems to be typical of the arcane nature of ggplot invocations. Did you
try tmap? Almost all the loaded&attached packages are superfluous too.

Roger

>
> The code I am using and the coordinates data is below. What am I doing
> wrong? Any help/suggestions will be highly appreciated.
>
> library(data.table)
> library(foreign)
> library(readstata13)
> library(rgdal)
> library(maptools)
> library(ggplot2)
> library(dplyr)
>
> load(url("
> https://github.com/valentinitnelav/RandomScripts/blob/master/NaturalEarth.RData?raw=true
> "))
>
> PROJ <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84
> +units=m +no_defs"
>
> NE_countries_rob  <- spTransform(NE_countries, CRSobj = PROJ)
> NE_graticules_rob <- spTransform(NE_graticules, CRSobj = PROJ)
> NE_box_rob        <- spTransform(NE_box, CRSobj = PROJ)
>
> # project long-lat coordinates for graticule label data frames
> # (two extra columns with projected XY are created)
> prj.coord <- project(cbind(lbl.Y$lon, lbl.Y$lat), proj=PROJ)
> lbl.Y.prj <- cbind(prj.coord, lbl.Y)
> names(lbl.Y.prj)[1:2] <- c("X.prj","Y.prj")
>
> prj.coord <- project(cbind(lbl.X$lon, lbl.X$lat), proj=PROJ)
> lbl.X.prj <- cbind(prj.coord, lbl.X)
> names(lbl.X.prj)[1:2] <- c("X.prj","Y.prj")
>
> m <- ggplot() +
>  # add Natural Earth countries projected to Robinson, give black border
> and fill with gray
>  geom_polygon(data=NE_countries_rob, aes(long,lat, group=group),
> colour="black", fill="gray80", size = 0.25) +
>  # Note: "Regions defined for each Polygons" warning has to do with
> fortify transformation. Might get deprecated in future!
>  # alternatively, use use map_data(NE_countries) to transform to data
> frame and then use project() to change to desired projection.
>  # add Natural Earth box projected to Robinson
>  geom_polygon(data=NE_box_rob, aes(x=long, y=lat), colour="black",
> fill="transparent", size = 0.25) +
>  # add graticules projected to Robinson
>  geom_path(data=NE_graticules_rob, aes(long, lat, group=group),
> linetype="dotted", color="grey50", size = 0.25) +
>  # add graticule labels - latitude and longitude
>  geom_text(data = lbl.Y.prj, aes(x = X.prj, y = Y.prj, label = lbl),
> color="grey50", size=2) +
>  geom_text(data = lbl.X.prj, aes(x = X.prj, y = Y.prj, label = lbl),
> color="grey50", size=2) +
>  # the default, ratio = 1 in coord_fixed ensures that one unit on the
> x-axis is the same length as one unit on the y-axis
>  coord_fixed(ratio = 1) +
>  geom_point(data=df,
>             aes(x=lon, y=lat), colour="Deep Pink",
>             fill="Pink",pch=21, size=2, alpha=I(0.7))
>  # remove the background and default gridlines
>  theme_void()
>
> ## coordinates dataframe
> dput(df)
> structure(list(lon = c(2.67569724621467, 17.5766416259819,
> 28.4126232192772,
> 23.8147674538232, 29.8917589327105), lat = c(28.1503115976162,
> -12.3388787380201, 9.78891068739477, -22.1873831176644, -3.36546931479253
> )), class = "data.frame", row.names = c(NA, -5L))
>
> Sincerely,
>
> Milu
>
> [[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]
https://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

Plot coordinates on world map with Robinson CRS - ggplot2

Mon, 02/18/2019 - 02:57
Dear all,

I am trying to plot coordinates on a world map with Robinson CRS. While the
world map is generated without any issues, when I try to plot the points -
I only get a single point.

The code I am using and the coordinates data is below. What am I doing
wrong? Any help/suggestions will be highly appreciated.

library(data.table)
library(foreign)
library(readstata13)
library(rgdal)
library(maptools)
library(ggplot2)
library(dplyr)

load(url("
https://github.com/valentinitnelav/RandomScripts/blob/master/NaturalEarth.RData?raw=true
"))

PROJ <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84
+units=m +no_defs"

NE_countries_rob  <- spTransform(NE_countries, CRSobj = PROJ)
NE_graticules_rob <- spTransform(NE_graticules, CRSobj = PROJ)
NE_box_rob        <- spTransform(NE_box, CRSobj = PROJ)

# project long-lat coordinates for graticule label data frames
# (two extra columns with projected XY are created)
prj.coord <- project(cbind(lbl.Y$lon, lbl.Y$lat), proj=PROJ)
lbl.Y.prj <- cbind(prj.coord, lbl.Y)
names(lbl.Y.prj)[1:2] <- c("X.prj","Y.prj")

prj.coord <- project(cbind(lbl.X$lon, lbl.X$lat), proj=PROJ)
lbl.X.prj <- cbind(prj.coord, lbl.X)
names(lbl.X.prj)[1:2] <- c("X.prj","Y.prj")

m <- ggplot() +
  # add Natural Earth countries projected to Robinson, give black border
and fill with gray
  geom_polygon(data=NE_countries_rob, aes(long,lat, group=group),
colour="black", fill="gray80", size = 0.25) +
  # Note: "Regions defined for each Polygons" warning has to do with
fortify transformation. Might get deprecated in future!
  # alternatively, use use map_data(NE_countries) to transform to data
frame and then use project() to change to desired projection.
  # add Natural Earth box projected to Robinson
  geom_polygon(data=NE_box_rob, aes(x=long, y=lat), colour="black",
fill="transparent", size = 0.25) +
  # add graticules projected to Robinson
  geom_path(data=NE_graticules_rob, aes(long, lat, group=group),
linetype="dotted", color="grey50", size = 0.25) +
  # add graticule labels - latitude and longitude
  geom_text(data = lbl.Y.prj, aes(x = X.prj, y = Y.prj, label = lbl),
color="grey50", size=2) +
  geom_text(data = lbl.X.prj, aes(x = X.prj, y = Y.prj, label = lbl),
color="grey50", size=2) +
  # the default, ratio = 1 in coord_fixed ensures that one unit on the
x-axis is the same length as one unit on the y-axis
  coord_fixed(ratio = 1) +
  geom_point(data=df,
             aes(x=lon, y=lat), colour="Deep Pink",
             fill="Pink",pch=21, size=2, alpha=I(0.7))
  # remove the background and default gridlines
  theme_void()

## coordinates dataframe
dput(df)
structure(list(lon = c(2.67569724621467, 17.5766416259819,
28.4126232192772,
23.8147674538232, 29.8917589327105), lat = c(28.1503115976162,
-12.3388787380201, 9.78891068739477, -22.1873831176644, -3.36546931479253
)), class = "data.frame", row.names = c(NA, -5L))

Sincerely,

Milu

        [[alternative HTML version deleted]]

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

Euclidean distance raster

Sun, 02/17/2019 - 20:46
Hi all,

I want to create a nearest distance raster in R. As far as I know, the
common practice is to compute the distance matrix of all the points in the
mask raster to all the features in the vector layer, then calculate the min
distance for each point. But this makes a huge (m by n) matrix and it
involves a lot of computation which is too slow.

To go around the computation, I wrote the following function that first
gets the index of the nearest feature to each point in the mask raster,
then only calculates the distance between each point and its nearest
feature. This solution is 3 times faster than the previous one, but still
not fast enough.

I wonder if anyone knows a better solution to this problem.

Thank you in advance,
Roozbeh Valavi

rasterDistance <- function(feature, rastermask){
  require(raster)
  require(sf)
  require(progress)
  p <- st_as_sf(rasterToPoints(rastermask, spatial = TRUE))
  p$indx <- st_nearest_feature(st_geometry(p), st_geometry(feature))
  pb <- progress::progress_bar$new(format = " Progress [:bar] :percent in
:elapsed",
                                   total=nrow(p), clear=FALSE, width=75) #
add progress bar
  for(i in 1:nrow(p)){
    p$dist[i] <- st_distance(st_geometry(p[i,]),
st_geometry(feature[p$indx[i],]))
    pb$tick() # update progress bar
  }
  output <- raster::rasterize(p, rastermask, field = "dist")
  return(output)
}

--
*Roozbeh Valavi*
PhD Candidate
The Quantitative & Applied Ecology Group <http://qaeco.com/>
School of BioSciences | Faculty of Science
The University of Melbourne, VIC 3010, Australia
Mobile: +61 423 283 238

        [[alternative HTML version deleted]]

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

Error with running bfastSpatial on MODIS data

Sun, 02/17/2019 - 06:35
Dear all,
I got error with running bfastSpatial on MODIS data for period of 2001-2019
. What could be my problem here?

###########

library(bfastSpatial)



##########

#MODIS data for period of 2001-2019

  n <- timeStackMODIS(list.files(path = dir, pattern = ".tif"))

   t <- system.time(

    bfm <- bfmSpatial(n, start=c(2009, 1), order=1)

 )



 # extract change raster

 change <- raster(bfm, 1)

months <- changeMonth(change)



###### Error

Error in setValues(r, values = methods::callGeneric(getValues(e1), e2)) :

  length(values) is not equal to ncell(x), or to 1

        [[alternative HTML version deleted]]

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

Re: Is there a way to merge different stacks?

Wed, 02/13/2019 - 09:28
Thank you for your answer.

Lara

Hugo Costa <[hidden email]> escreveu no dia quarta, 13/02/2019 à(s)
14:21:

> You are possibly looking for function mosaic (raster package).
> Hugo
>
> Lara Silva <[hidden email]> escreveu no dia quarta, 13/02/2019
> à(s) 14:47:
>
>> Hello everybody,
>> I am having a problem to merge several stacks. My stacks have different
>> extensions ... It is a problem.
>>
>> If possible, I want to merge several (5 stacks) asc grids into in a
>> "global
>> grid/stack"
>> How can I do it? "Merging", "joining" ...?
>> Any sugestion?
>>
>> Thanks,
>>
>> Lara Silva
>>
>> My outputs:
>>
>> 1st stack
>>
>> > grids_climatic_SMP <- list.files(pattern='asc', full.names=T)# All "asc"
>> files in folder
>> > climatic_SMP<- stack(grids_climatic_SMP)# Make a "stack" of EGV #Install
>> "stack"
>> > climatic_SMP
>> class       : RasterStack
>> dimensions  : 901, 901, 811801, 20  (nrow, ncol, ncell, nlayers)
>> resolution  : 100, 100  (x, y)
>> extent      : 586974, 677074, 4140847, 4230947  (xmin, xmax, ymin, ymax)
>> coord. ref. : NA
>> names       : bio1_SMP, bio10_SMP, bio11_SMP, bio12_SMP, bio13_SMP,
>> bio14_SMP, bio15_SMP, bio16_SMP, bio17_SMP, bio18_SMP, bio19_SMP,
>> bio2_SMP,
>> bio3_SMP, bio4_SMP, bio4a_SMP, ...
>>
>> 2nd stack
>>
>> > grids_climatic_PP <- list.files(pattern='asc', full.names=T)# All "asc"
>> files in folder
>> > climatic_PP<- stack(grids_climatic_PP)# Make a "stack" of EGV #Install
>> "stack"
>> > climatic_PP
>> class       : RasterStack
>> dimensions  : 601, 601, 361201, 20  (nrow, ncol, ncell, nlayers)
>> resolution  : 100, 100  (x, y)
>> extent      : 358603, 418703, 4229376, 4289476  (xmin, xmax, ymin, ymax)
>> coord. ref. : NA
>> names       : bio1_PP, bio10_PP, bio11_PP, bio12_PP, bio13_PP, bio14_PP,
>> bio15_PP, bio16_PP, bio17_PP, bio18_PP, bio19_PP, bio2_PP, bio3_PP,
>> bio4_PP, bio4a_PP, ...
>>
>> 3rd stack
>>
>> > grids_climatic_TP <- list.files(pattern='asc', full.names=T)# All "asc"
>> files in folder
>> > climatic_TP<- stack(grids_climatic_TP)# Make a "stack" of EGV #Install
>> "stack"
>> > climatic_TP
>> class       : RasterStack
>> dimensions  : 277, 412, 114124, 20  (nrow, ncol, ncell, nlayers)
>> resolution  : 100, 100  (x, y)
>> extent      : 461214, 502414, 4272892, 4300592  (xmin, xmax, ymin, ymax)
>> coord. ref. : NA
>> names       : bio1_TP, bio10_TP, bio11_TP, bio12_TP, bio13_TP, bio14_TP,
>> bio15_TP, bio16_TP, bio17_TP, bio18_TP, bio19_TP, bio2_TP, bio3_TP,
>> bio4_TP, bio4a_TP, ...
>>
>>
>> 4th stack
>>
>> > grids_climatic_FP <- list.files(pattern='asc', full.names=T)# All "asc"
>> files in folder
>> > climatic_FP<- stack(grids_climatic_FP)# Make a "stack" of EGV #Install
>> "stack"
>> > climatic_FP
>> class       : RasterStack
>> dimensions  : 301, 301, 90601, 20  (nrow, ncol, ncell, nlayers)
>> resolution  : 100, 100  (x, y)
>> extent      : 336088, 366188, 4256440, 4286540  (xmin, xmax, ymin, ymax)
>> coord. ref. : NA
>> names       : bio1_FP, bio10_FP, bio11_FP, bio12_FP, bio13_FP, bio14_FP,
>> bio15_FP, bio16_FP, bio17_FP, bio18_FP, bio19_FP, bio2_FP, bio3_FP,
>> bio4_FP, bio4a_FP, ...
>>
>>
>> 5th stack
>>
>> > grids_climatic_SJP <- list.files(pattern='asc', full.names=T)# All "asc"
>> files in folder
>> > climatic_SJP<- stack(grids_climatic_SJP)# Make a "stack" of EGV #Install
>> "stack"
>> > climatic_SJP
>> class       : RasterStack
>> dimensions  : 602, 602, 362404, 20  (nrow, ncol, ncell, nlayers)
>> resolution  : 100, 100  (x, y)
>> extent      : 381100, 441300, 4249487, 4309687  (xmin, xmax, ymin, ymax)
>> coord. ref. : NA
>> names       : bio1_SJP, bio10_SJP, bio11_SJP, bio12_SJP, bio13_SJP,
>> bio14_SJP, bio15_SJP, bio16_SJP, bio17_SJP, bio18_SJP, bio19_SJP,
>> bio2_SJP,
>> bio3_SJP, bio4_SJP, bio4a_SJP, ...
>>
>> >
>>
>>         [[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: Is there a way to merge different stacks?

Wed, 02/13/2019 - 09:21
You are possibly looking for function mosaic (raster package).
Hugo

Lara Silva <[hidden email]> escreveu no dia quarta, 13/02/2019
à(s) 14:47:

> Hello everybody,
> I am having a problem to merge several stacks. My stacks have different
> extensions ... It is a problem.
>
> If possible, I want to merge several (5 stacks) asc grids into in a "global
> grid/stack"
> How can I do it? "Merging", "joining" ...?
> Any sugestion?
>
> Thanks,
>
> Lara Silva
>
> My outputs:
>
> 1st stack
>
> > grids_climatic_SMP <- list.files(pattern='asc', full.names=T)# All "asc"
> files in folder
> > climatic_SMP<- stack(grids_climatic_SMP)# Make a "stack" of EGV #Install
> "stack"
> > climatic_SMP
> class       : RasterStack
> dimensions  : 901, 901, 811801, 20  (nrow, ncol, ncell, nlayers)
> resolution  : 100, 100  (x, y)
> extent      : 586974, 677074, 4140847, 4230947  (xmin, xmax, ymin, ymax)
> coord. ref. : NA
> names       : bio1_SMP, bio10_SMP, bio11_SMP, bio12_SMP, bio13_SMP,
> bio14_SMP, bio15_SMP, bio16_SMP, bio17_SMP, bio18_SMP, bio19_SMP, bio2_SMP,
> bio3_SMP, bio4_SMP, bio4a_SMP, ...
>
> 2nd stack
>
> > grids_climatic_PP <- list.files(pattern='asc', full.names=T)# All "asc"
> files in folder
> > climatic_PP<- stack(grids_climatic_PP)# Make a "stack" of EGV #Install
> "stack"
> > climatic_PP
> class       : RasterStack
> dimensions  : 601, 601, 361201, 20  (nrow, ncol, ncell, nlayers)
> resolution  : 100, 100  (x, y)
> extent      : 358603, 418703, 4229376, 4289476  (xmin, xmax, ymin, ymax)
> coord. ref. : NA
> names       : bio1_PP, bio10_PP, bio11_PP, bio12_PP, bio13_PP, bio14_PP,
> bio15_PP, bio16_PP, bio17_PP, bio18_PP, bio19_PP, bio2_PP, bio3_PP,
> bio4_PP, bio4a_PP, ...
>
> 3rd stack
>
> > grids_climatic_TP <- list.files(pattern='asc', full.names=T)# All "asc"
> files in folder
> > climatic_TP<- stack(grids_climatic_TP)# Make a "stack" of EGV #Install
> "stack"
> > climatic_TP
> class       : RasterStack
> dimensions  : 277, 412, 114124, 20  (nrow, ncol, ncell, nlayers)
> resolution  : 100, 100  (x, y)
> extent      : 461214, 502414, 4272892, 4300592  (xmin, xmax, ymin, ymax)
> coord. ref. : NA
> names       : bio1_TP, bio10_TP, bio11_TP, bio12_TP, bio13_TP, bio14_TP,
> bio15_TP, bio16_TP, bio17_TP, bio18_TP, bio19_TP, bio2_TP, bio3_TP,
> bio4_TP, bio4a_TP, ...
>
>
> 4th stack
>
> > grids_climatic_FP <- list.files(pattern='asc', full.names=T)# All "asc"
> files in folder
> > climatic_FP<- stack(grids_climatic_FP)# Make a "stack" of EGV #Install
> "stack"
> > climatic_FP
> class       : RasterStack
> dimensions  : 301, 301, 90601, 20  (nrow, ncol, ncell, nlayers)
> resolution  : 100, 100  (x, y)
> extent      : 336088, 366188, 4256440, 4286540  (xmin, xmax, ymin, ymax)
> coord. ref. : NA
> names       : bio1_FP, bio10_FP, bio11_FP, bio12_FP, bio13_FP, bio14_FP,
> bio15_FP, bio16_FP, bio17_FP, bio18_FP, bio19_FP, bio2_FP, bio3_FP,
> bio4_FP, bio4a_FP, ...
>
>
> 5th stack
>
> > grids_climatic_SJP <- list.files(pattern='asc', full.names=T)# All "asc"
> files in folder
> > climatic_SJP<- stack(grids_climatic_SJP)# Make a "stack" of EGV #Install
> "stack"
> > climatic_SJP
> class       : RasterStack
> dimensions  : 602, 602, 362404, 20  (nrow, ncol, ncell, nlayers)
> resolution  : 100, 100  (x, y)
> extent      : 381100, 441300, 4249487, 4309687  (xmin, xmax, ymin, ymax)
> coord. ref. : NA
> names       : bio1_SJP, bio10_SJP, bio11_SJP, bio12_SJP, bio13_SJP,
> bio14_SJP, bio15_SJP, bio16_SJP, bio17_SJP, bio18_SJP, bio19_SJP, bio2_SJP,
> bio3_SJP, bio4_SJP, bio4a_SJP, ...
>
> >
>
>         [[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

Is there a way to merge different stacks?

Wed, 02/13/2019 - 08:47
Hello everybody,
I am having a problem to merge several stacks. My stacks have different
extensions ... It is a problem.

If possible, I want to merge several (5 stacks) asc grids into in a "global
grid/stack"
How can I do it? "Merging", "joining" ...?
Any sugestion?

Thanks,

Lara Silva

My outputs:

1st stack

> grids_climatic_SMP <- list.files(pattern='asc', full.names=T)# All "asc"
files in folder
> climatic_SMP<- stack(grids_climatic_SMP)# Make a "stack" of EGV #Install
"stack"
> climatic_SMP
class       : RasterStack
dimensions  : 901, 901, 811801, 20  (nrow, ncol, ncell, nlayers)
resolution  : 100, 100  (x, y)
extent      : 586974, 677074, 4140847, 4230947  (xmin, xmax, ymin, ymax)
coord. ref. : NA
names       : bio1_SMP, bio10_SMP, bio11_SMP, bio12_SMP, bio13_SMP,
bio14_SMP, bio15_SMP, bio16_SMP, bio17_SMP, bio18_SMP, bio19_SMP, bio2_SMP,
bio3_SMP, bio4_SMP, bio4a_SMP, ...

2nd stack

> grids_climatic_PP <- list.files(pattern='asc', full.names=T)# All "asc"
files in folder
> climatic_PP<- stack(grids_climatic_PP)# Make a "stack" of EGV #Install
"stack"
> climatic_PP
class       : RasterStack
dimensions  : 601, 601, 361201, 20  (nrow, ncol, ncell, nlayers)
resolution  : 100, 100  (x, y)
extent      : 358603, 418703, 4229376, 4289476  (xmin, xmax, ymin, ymax)
coord. ref. : NA
names       : bio1_PP, bio10_PP, bio11_PP, bio12_PP, bio13_PP, bio14_PP,
bio15_PP, bio16_PP, bio17_PP, bio18_PP, bio19_PP, bio2_PP, bio3_PP,
bio4_PP, bio4a_PP, ...

3rd stack

> grids_climatic_TP <- list.files(pattern='asc', full.names=T)# All "asc"
files in folder
> climatic_TP<- stack(grids_climatic_TP)# Make a "stack" of EGV #Install
"stack"
> climatic_TP
class       : RasterStack
dimensions  : 277, 412, 114124, 20  (nrow, ncol, ncell, nlayers)
resolution  : 100, 100  (x, y)
extent      : 461214, 502414, 4272892, 4300592  (xmin, xmax, ymin, ymax)
coord. ref. : NA
names       : bio1_TP, bio10_TP, bio11_TP, bio12_TP, bio13_TP, bio14_TP,
bio15_TP, bio16_TP, bio17_TP, bio18_TP, bio19_TP, bio2_TP, bio3_TP,
bio4_TP, bio4a_TP, ...


4th stack

> grids_climatic_FP <- list.files(pattern='asc', full.names=T)# All "asc"
files in folder
> climatic_FP<- stack(grids_climatic_FP)# Make a "stack" of EGV #Install
"stack"
> climatic_FP
class       : RasterStack
dimensions  : 301, 301, 90601, 20  (nrow, ncol, ncell, nlayers)
resolution  : 100, 100  (x, y)
extent      : 336088, 366188, 4256440, 4286540  (xmin, xmax, ymin, ymax)
coord. ref. : NA
names       : bio1_FP, bio10_FP, bio11_FP, bio12_FP, bio13_FP, bio14_FP,
bio15_FP, bio16_FP, bio17_FP, bio18_FP, bio19_FP, bio2_FP, bio3_FP,
bio4_FP, bio4a_FP, ...


5th stack

> grids_climatic_SJP <- list.files(pattern='asc', full.names=T)# All "asc"
files in folder
> climatic_SJP<- stack(grids_climatic_SJP)# Make a "stack" of EGV #Install
"stack"
> climatic_SJP
class       : RasterStack
dimensions  : 602, 602, 362404, 20  (nrow, ncol, ncell, nlayers)
resolution  : 100, 100  (x, y)
extent      : 381100, 441300, 4249487, 4309687  (xmin, xmax, ymin, ymax)
coord. ref. : NA
names       : bio1_SJP, bio10_SJP, bio11_SJP, bio12_SJP, bio13_SJP,
bio14_SJP, bio15_SJP, bio16_SJP, bio17_SJP, bio18_SJP, bio19_SJP, bio2_SJP,
bio3_SJP, bio4_SJP, bio4a_SJP, ...

>

        [[alternative HTML version deleted]]

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

Pages