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

Getis Ord G test in 3 dimensions

Wed, 03/15/2017 - 11:43
> Dear all

> I am working on a 3D home range estimator, we want to calculated home ranges and core areas on primates that were followed sleeping site to sleeping site. We collected GPS locations and animals height at equal intervals along the day.
> We already calculates home ranges and core area in 2D using Delaunay triangulation and then applying a hot spot analysis in ArcGis( Getis Ord G*) using the perimeter of the triangles. Then, to calculate it on 3D we first created a 3D convex hull formed by tetrahedrons using animals locations (coordinates in meters)  and height.

> we are wondering if it is possible to use the function Local G or Global G test (package Spdep) for 3D locatios (X,Y,Z) as we previously  did in 2D (Hotspot with rendering, ArgGis)

Thanks in advance
>
>
> Juanma
                       
 
Juan Manuel José Domínguez
Postdoctoral Researcher        
                                                 
Conservation Ecology Program
School of Bioresources & Technology
King Mongkut´s University of Technology Thonburi

49 Soi Tienthalay 25
Bangkhuntien-Chaithalay Road
Thakham, Bangkhuntien
Bangkok 10150
Thailand
http://cons-ecol-kmutt.weebly.com/
        [[alternative HTML version deleted]]

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

Re: loop memory problem

Wed, 03/15/2017 - 11:37
Hi again,

I read in the guide of foreach that include the packages was a good idea.
Now the parallel loop is runing but at the end I got only one spatial line
(the first) I don't understand why is not adding the other spatiallines.

Any idea?

###3)dopar%
library(doSNOW)
registerDoSEQ()
registerDoParallel(cores=10)
getDoParWorkers()
line2<- gdistance::shortestPath(trCostS16, pos@coords[13,],pos@coords[12,],
output="SpatialLines")#

x<-foreach(i=2:13,.packages=c("gdistance","doParallel","foreach","base","sp"))
%dopar% {
  nextSegment2=gdistance::shortestPath(trCostS16,
pos@coords[i,],pos@coords[i+1,],
output="SpatialLines")
  line2 =nextSegment2+line2
  print(i)}

stopCluster(cl)
registerDoSEQ()

Thanks in advance,
Marta
####files
https://drive.google.com/open?id=0BwqSBe1Yq-FBUWVBOUdvaThEU1k

        [[alternative HTML version deleted]]

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

Re: Same mask for two rasters, but obtaining different extents in R

Wed, 03/15/2017 - 11:21
I re-checked the original extent and here it is: This is data obtained from
WorldClim.

> bio5
class       : RasterLayer
dimensions  : 3600, 8640, 31104000  (nrow, ncol, ncell)
resolution  : 0.04166667, 0.04166667  (x, y)
extent      : -180, 180, -60, *90.00001 * (xmin, xmax, ymin, ymax) #This
varies at the 6th decimal place.
coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs
data source : C:\Users\rameshv\Downloads\Climate
Stability\Data_LGM_Present\Present\1_RawData\bio_5
names       : bio_5
values      : -86, 489  (min, max)
attributes  :
        ID COUNT
 from: -86     1
 to  : 489    38


> lg5
class       : RasterLayer
dimensions  : 3600, 8640, 31104000  (nrow, ncol, ncell)
resolution  : 0.04166667, 0.04166667  (x, y)
extent      : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs
data source : C:\Users\rameshv\Downloads\Climate
Stability\Data_LGM_Present\LGM\1_RawData\cclgmbi5.tif
names       : cclgmbi5
values      : -204, 519  (min, max)

What do you suggest?

Also, any suggestions on aligning them before I crop them? I wouldn't want
to projectRaster using one of the rasters, as that creates extra number of
rows (when I coerce the raster to a Spatial Pixels Data Frame).



On Wed, Mar 15, 2017 at 1:13 PM, Sarah Goslee <[hidden email]>
wrote:

> The most likely explanation seems to me that the rasters were not
> aligned before cropping, so they are not aligned after cropping.
>
> Have you checked the original extent and resolution?
>
> Sarah
>
> On Tue, Mar 14, 2017 at 5:31 PM, Vijay Ramesh via R-sig-Geo
> <[hidden email]> wrote:
> > I loaded a shapefile in R, and cropped and masked it with the same mask,
> > but I get different extents. Any suggestions?
> >
> > Code below:
> >
> >     library(raster)
> >     library(rgdal)
> >     library(GISTools)
> >     library(sp)
> >     library(maptools)
> >
> >     ##Loading the first mask
> >     Mask <- readOGR("C:\\Users\\rameshv\\Downloads\\Climate
> > Stability\\SPV_Mask\\mask.shp")
> >     proj4string(Mask) <- CRS("+init=epsg:4326")
> >
> >     #Loading the Rasters For Present and LGM and clipping them to the
> right
> > extent
> >     bio5 <- raster("C:\\Users\\rameshv\\Downloads\\Climate
> > Stability\\Data_LGM_Present\\Present\\1_RawData\\bio_5")
> >     proj4string(bio5) <- CRS("+init=epsg:4326")
> >     lg5 <- raster("C:\\Users\\rameshv\\Downloads\\Climate
> > Stability\\Data_LGM_Present\\LGM\\1_RawData\\cclgmbi5.tif")
> >     proj4string(lg5) <- CRS("+init=epsg:4326")
> >
> >     ##Cropping by using the Crop and Mask functions
> >     cr<- crop(bio5, extent(Mask))
> >     bio5 <- mask(cr, Mask)
> >
> >     cr2<- crop(lg5, extent(Mask))
> >     lg5<- mask(cr2, Mask)
> >
> >     > bio5
> >      class       : RasterLayer
> >      dimensions  : 339, 246, 83394  (nrow, ncol, ncell)
> >      resolution  : 0.04166667, 0.04166667  (x, y)
> >      extent      : 72.45835, 82.70835, 8.083337, 22.20834  (xmin, xmax,
> > ymin, ymax)
> >      coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs
> > +ellps=WGS84 +towgs84=0,0,0
> >      data source : in memory
> >      names       : bio_5
> >      values      : 207, 432  (min, max)
> >      attributes  :
> >         ID COUNT
> >      from: -86     1
> >      to  : 489    38
> >
> >      > lg5
> >      class       : RasterLayer
> >      dimensions  : 339, 246, 83394  (nrow, ncol, ncell)
> >      resolution  : 0.04166667, 0.04166667  (x, y)
> >      extent      : 72.45833, 82.70833, 8.083333, 22.20833  (xmin, xmax,
> > ymin, ymax)
> >      coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs
> > +ellps=WGS84 +towgs84=0,0,0
> >      data source : in memory
> >      names       : cclgmbi5
> >      values      : 187, 392  (min, max)
> >
> >
>
> --
> Sarah Goslee
> http://www.functionaldiversity.org
>


--
Data Manager,
Barbara Han Lab,
Cary Institute of Ecosystem Studies,
2801 Sharon Turnpike, Millbrook, NY 12545
Lab Site : http://www.hanlab.science/
Personal Site : http://evolecol.weebly.com/
Phone : (845)-677-7600 Ext: 241

        [[alternative HTML version deleted]]

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

Re: Same mask for two rasters, but obtaining different extents in R

Wed, 03/15/2017 - 11:13
The most likely explanation seems to me that the rasters were not
aligned before cropping, so they are not aligned after cropping.

Have you checked the original extent and resolution?

Sarah

On Tue, Mar 14, 2017 at 5:31 PM, Vijay Ramesh via R-sig-Geo
<[hidden email]> wrote:
> I loaded a shapefile in R, and cropped and masked it with the same mask,
> but I get different extents. Any suggestions?
>
> Code below:
>
>     library(raster)
>     library(rgdal)
>     library(GISTools)
>     library(sp)
>     library(maptools)
>
>     ##Loading the first mask
>     Mask <- readOGR("C:\\Users\\rameshv\\Downloads\\Climate
> Stability\\SPV_Mask\\mask.shp")
>     proj4string(Mask) <- CRS("+init=epsg:4326")
>
>     #Loading the Rasters For Present and LGM and clipping them to the right
> extent
>     bio5 <- raster("C:\\Users\\rameshv\\Downloads\\Climate
> Stability\\Data_LGM_Present\\Present\\1_RawData\\bio_5")
>     proj4string(bio5) <- CRS("+init=epsg:4326")
>     lg5 <- raster("C:\\Users\\rameshv\\Downloads\\Climate
> Stability\\Data_LGM_Present\\LGM\\1_RawData\\cclgmbi5.tif")
>     proj4string(lg5) <- CRS("+init=epsg:4326")
>
>     ##Cropping by using the Crop and Mask functions
>     cr<- crop(bio5, extent(Mask))
>     bio5 <- mask(cr, Mask)
>
>     cr2<- crop(lg5, extent(Mask))
>     lg5<- mask(cr2, Mask)
>
>     > bio5
>      class       : RasterLayer
>      dimensions  : 339, 246, 83394  (nrow, ncol, ncell)
>      resolution  : 0.04166667, 0.04166667  (x, y)
>      extent      : 72.45835, 82.70835, 8.083337, 22.20834  (xmin, xmax,
> ymin, ymax)
>      coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs
> +ellps=WGS84 +towgs84=0,0,0
>      data source : in memory
>      names       : bio_5
>      values      : 207, 432  (min, max)
>      attributes  :
>         ID COUNT
>      from: -86     1
>      to  : 489    38
>
>      > lg5
>      class       : RasterLayer
>      dimensions  : 339, 246, 83394  (nrow, ncol, ncell)
>      resolution  : 0.04166667, 0.04166667  (x, y)
>      extent      : 72.45833, 82.70833, 8.083333, 22.20833  (xmin, xmax,
> ymin, ymax)
>      coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs
> +ellps=WGS84 +towgs84=0,0,0
>      data source : in memory
>      names       : cclgmbi5
>      values      : 187, 392  (min, max)
>
>
--
Sarah Goslee
http://www.functionaldiversity.org

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

Re: Same mask for two rasters, but obtaining different extents in R

Wed, 03/15/2017 - 07:58
Hi,

Might you show us what bio5 and lg5 look like before you do the masking?

Cheers,
Ben


> On Mar 14, 2017, at 5:31 PM, Vijay Ramesh via R-sig-Geo <[hidden email]> wrote:
>
> I loaded a shapefile in R, and cropped and masked it with the same mask,
> but I get different extents. Any suggestions?
>
> Code below:
>
>    library(raster)
>    library(rgdal)
>    library(GISTools)
>    library(sp)
>    library(maptools)
>
>    ##Loading the first mask
>    Mask <- readOGR("C:\\Users\\rameshv\\Downloads\\Climate
> Stability\\SPV_Mask\\mask.shp")
>    proj4string(Mask) <- CRS("+init=epsg:4326")
>
>    #Loading the Rasters For Present and LGM and clipping them to the right
> extent
>    bio5 <- raster("C:\\Users\\rameshv\\Downloads\\Climate
> Stability\\Data_LGM_Present\\Present\\1_RawData\\bio_5")
>    proj4string(bio5) <- CRS("+init=epsg:4326")
>    lg5 <- raster("C:\\Users\\rameshv\\Downloads\\Climate
> Stability\\Data_LGM_Present\\LGM\\1_RawData\\cclgmbi5.tif")
>    proj4string(lg5) <- CRS("+init=epsg:4326")
>
>    ##Cropping by using the Crop and Mask functions
>    cr<- crop(bio5, extent(Mask))
>    bio5 <- mask(cr, Mask)
>
>    cr2<- crop(lg5, extent(Mask))
>    lg5<- mask(cr2, Mask)
>
>> bio5
>     class       : RasterLayer
>     dimensions  : 339, 246, 83394  (nrow, ncol, ncell)
>     resolution  : 0.04166667, 0.04166667  (x, y)
>     extent      : 72.45835, 82.70835, 8.083337, 22.20834  (xmin, xmax,
> ymin, ymax)
>     coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs
> +ellps=WGS84 +towgs84=0,0,0
>     data source : in memory
>     names       : bio_5
>     values      : 207, 432  (min, max)
>     attributes  :
>        ID COUNT
>     from: -86     1
>     to  : 489    38
>
>> lg5
>     class       : RasterLayer
>     dimensions  : 339, 246, 83394  (nrow, ncol, ncell)
>     resolution  : 0.04166667, 0.04166667  (x, y)
>     extent      : 72.45833, 82.70833, 8.083333, 22.20833  (xmin, xmax,
> ymin, ymax)
>     coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs
> +ellps=WGS84 +towgs84=0,0,0
>     data source : in memory
>     names       : cclgmbi5
>     values      : 187, 392  (min, max)
>
>
> --
> Data Manager,
> Barbara Han Lab,
> Cary Institute of Ecosystem Studies,
> 2801 Sharon Turnpike, Millbrook, NY 12545
> Lab Site : http://www.hanlab.science/
> Personal Site : http://evolecol.weebly.com/
> Phone : (845)-677-7600 Ext: 241
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Ben Tupper
Bigelow Laboratory for Ocean Sciences
60 Bigelow Drive, P.O. Box 380
East Boothbay, Maine 04544
http://www.bigelow.org

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

Same mask for two rasters, but obtaining different extents in R

Tue, 03/14/2017 - 15:31
I loaded a shapefile in R, and cropped and masked it with the same mask,
but I get different extents. Any suggestions?

Code below:

    library(raster)
    library(rgdal)
    library(GISTools)
    library(sp)
    library(maptools)

    ##Loading the first mask
    Mask <- readOGR("C:\\Users\\rameshv\\Downloads\\Climate
Stability\\SPV_Mask\\mask.shp")
    proj4string(Mask) <- CRS("+init=epsg:4326")

    #Loading the Rasters For Present and LGM and clipping them to the right
extent
    bio5 <- raster("C:\\Users\\rameshv\\Downloads\\Climate
Stability\\Data_LGM_Present\\Present\\1_RawData\\bio_5")
    proj4string(bio5) <- CRS("+init=epsg:4326")
    lg5 <- raster("C:\\Users\\rameshv\\Downloads\\Climate
Stability\\Data_LGM_Present\\LGM\\1_RawData\\cclgmbi5.tif")
    proj4string(lg5) <- CRS("+init=epsg:4326")

    ##Cropping by using the Crop and Mask functions
    cr<- crop(bio5, extent(Mask))
    bio5 <- mask(cr, Mask)

    cr2<- crop(lg5, extent(Mask))
    lg5<- mask(cr2, Mask)

    > bio5
     class       : RasterLayer
     dimensions  : 339, 246, 83394  (nrow, ncol, ncell)
     resolution  : 0.04166667, 0.04166667  (x, y)
     extent      : 72.45835, 82.70835, 8.083337, 22.20834  (xmin, xmax,
ymin, ymax)
     coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs
+ellps=WGS84 +towgs84=0,0,0
     data source : in memory
     names       : bio_5
     values      : 207, 432  (min, max)
     attributes  :
        ID COUNT
     from: -86     1
     to  : 489    38

     > lg5
     class       : RasterLayer
     dimensions  : 339, 246, 83394  (nrow, ncol, ncell)
     resolution  : 0.04166667, 0.04166667  (x, y)
     extent      : 72.45833, 82.70833, 8.083333, 22.20833  (xmin, xmax,
ymin, ymax)
     coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs
+ellps=WGS84 +towgs84=0,0,0
     data source : in memory
     names       : cclgmbi5
     values      : 187, 392  (min, max)


--
Data Manager,
Barbara Han Lab,
Cary Institute of Ecosystem Studies,
2801 Sharon Turnpike, Millbrook, NY 12545
Lab Site : http://www.hanlab.science/
Personal Site : http://evolecol.weebly.com/
Phone : (845)-677-7600 Ext: 241

        [[alternative HTML version deleted]]

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

Re: loop memory problem

Tue, 03/14/2017 - 11:02
> #answer Marcel in r-sig-geo forum:
>> #1) maybe you can first make your calculations and afterwards plot the
>> result.
>> #It might be faster this way.
>>
>> #2) maybe you can use the parallel for each loop in order to use the
>> whole performance of your cpu
>>
>> #3) loops in R are really slow in general.. you could also think about
>> some fancy stuff like some R compiler or Rcpp package
>>
>> #4) in order to analyze your memory consumption, you should have a look
>> on your system resources (depending on your OS: task manager (win) or htop
>> (linux))
>>
>> #############
>> #1)MARCEL: maybe you can first make your calculations and afterwards plot
>> the result.
>> #It might be faster this way.
>> ##MARTA: I did that, but the difference between the calculations inside
>> the loop of afterwards is less than a second.
>> #
>>
>       #############
> #2) MARCEL:maybe you can use the parallel or each loop in order to use the
>> whole performance of your cpu
>> #MARTA:I  developed 5 different loops, the old one (2.A) it works
>> properly but is too slow. I followed your suggestion with the parallel
>> loops to increase the speed. I wrote a loop (2.B) with %do%, which is works
>> but the time is the #same as the old loop. The other loops with %doPar%,
>> didn't work at all. Only the simple (2.D) runs, but is not what I want. The
>> other two ( 2.C and 2.E) didn't run, they have errors. The 2.C "task 1
>> failed - "non-numeric argument #to binary operator"  and the 2.E "task 1
>> failed - "no method for coercing this S4 class to a vector".
>> #
>> #Any idea of how repair this loop?
>>
>> #########
>> #library
>> ########
>> library(parallel);library(foreach);library(doParallel);
>> library(gdistance);library(raster)
>> library(rgdal);library(rgeos)
>> library(sp)
>> #######
>> #data
>> ###########
>> path_data<-"E:/Q11/"
>> # cost raster sea and island
>> costa6Azo <- raster("E:/Q11/costa6Azo_projected.tif")
>> transitioncosta6Azo <- transition(costa6Azo, min, directions=16)#porque
>> min????
>> trCostS16 <- gdistance::geoCorrection(transitioncosta6Azo, type="c")
>> #points
>> boat <- read.table(paste0(path_data,"boat2905.csv"), header=TRUE,
>> sep=";", na.strings="NA", dec=".", strip.white=TRUE)
>> pos<-as.data.frame(cbind(boat$Lat1,boat$Long1,boat$Ref))
>> str(pos)
>> sp::coordinates(pos) <- ~V2+V1
>> sp::proj4string(pos) <-CRS("+proj=longlat +datum=WGS84 +no_defs")
>>
>> pos<-sp::spTransform(pos, CRS("+proj=utm +zone=26 +ellps=intl
>> +towgs84=-104,167,-38,0,0,0,0 +units=m +no_defs"))
>> # The aim of the loop: calculating the track of the whole sailing
>> path#############################################################
>> #
>> line<- gdistance::shortestPath(trCostS16, pos@coords[1,],pos@coords[2,],
>> output="SpatialLines")#
>> glength(1-2=4887.737m);(2-3=12590.01);(12-11=11360.39m);(12-13=9453.001m)
>> lines(line,col=5)
>> #2. A) old
>> loop##########################################################################################################################
>> #it works, but it's too slow
>>
>
> ## here we start with the for-loop
>> for (i in (seq(2,length(pos) - 1))) {
>>   # calculation of the rest of the segements
>>   nextSegment<- gdistance::shortestPath(trCostS16, pos@coords
>> [i,],pos@coords[i+1,], output="SpatialLines")
>>   # simple addition combines the single spatialline segements
>>   line <- nextSegment + line
>>   # we plot each new segment
>>   lines(nextSegment)
>> }
>> # note that we have now ten combined line features in this SpatialLines
>> object
>> line
>> gLength(line)#110747.2
>>
>> #2.B### new parallel
>> loop##################################################################################################3
>> #it works, but it's too slow##%do%
>> #
>> line<- gdistance::shortestPath(trCostS16, pos@coords[1,],pos@coords[2,],
>> output="SpatialLines")#
>> glength(1-2=4887.737m);(2-3=12590.01);(12-11=11360.39m);(12-13=9453.001m)
>>
>> x <- foreach(i=2:13) %do%
>>   {nextSegment<-gdistance::shortestPath(trCostS16, pos@coords
>> [i,],pos@coords[i-1,], output="SpatialLines")
>>     line <- nextSegment + line
>>  }
>> x      ##  it works!!!! 1+ 12 spatialLines!!
>> ##
>
> #     #2.C#parallel loop ##%dopar%

> #without success
>>
> registerDoParallel()
>
> getDoParWorkers()
>> line<- gdistance::shortestPath(trCostS16, pos@coords[1,],pos@coords[2,],
>> output="SpatialLines")#
>> glength(1-2=4887.737m);(2-3=12590.01);(12-11=11360.39m);(12-13=9453.001m)
>> #function
>> #
>> funMTM<-function(){
>> nextSegment<-gdistance::shortestPath(trCostS16, pos@coords[i,],pos@coords[i-1,],
>> output="SpatialLines")
>> line <- nextSegment + line
>> }
>> getDoParWorkers()
>> ptime <- system.time({
>>   result <- foreach(i=2:13) %dopar% funMTM()
>>   })
>> ptime#Error in funMTM() :
>> #task 1 failed - "non-numeric argument to binary operator"
>> #
>>
> #
> #2.D#loop %dopar% simple
>>
>        # it works, but I need the spatiallines output, not a list.

> #If I run the %dopar% only with the functions shortestPath without
>> increase the lines into the SpatialLines. I get a list with 13 features(
>> indidivual spatialLines). However, I need an SpatialLines object, with 13
>> spatiallines inside.
>> registerDoParallel()
>> registerDoSEQ()
>> registerDoParallel(cores=10)
>> getDoParWorkers()
>> system.time(foreach(i=2:13) %dopar% gdistance::shortestPath(trCostS16,
>> pos@coords[i,],pos@coords[i-1,], output="SpatialLines"))
>> #user  system elapsed
>> #0.74    0.69   54.81
>>
> #    #2.E#loop %dopar% with the function defined in the loop
      #without success:

> system.time(foreach(i=2:13) %dopar% {
>>   nextSegment=gdistance::shortestPath(trCostS16, pos@coords
>> [i,],pos@coords[i-1,], output="SpatialLines")
>>             line =nextSegment + line})
>> #Error in { :
>> #task 1 failed - "no method for coercing this S4 class to a vector"
>> stopCluster(cl)
>>
>> ##############
>> #3)MARCEL: loops in R are really slow in general.. you could also think
>> about some fancy stuff like some R compiler or Rcpp package
>> #MARTA##### Functions #####
>> # byte code compilation
>> library(compiler)
>> myfunc<-gdistance::shortestPath(trCostS16, pos@coords[1,],pos@coords[2,],
>> output="SpatialLines")
>> myFuncCmp <- cmpfun(myfunc)
>> system.time({
>>   output <- SpatialLines(LinesList = , 1, FUN=myFuncCmp)
>> })
>> #############
>>
>>       #############
> #4) MARCEL: in order to analyze your memory consumption, you should have a
>> look on your system resources (depending on your OS: task manager (win) or
>> htop (linux))
>> #MARTA:I run the loop with the task manager's window open, and never over
>> pass the 30% of the CPU's memory .
>> ####files
>> https://drive.google.com/open?id=0BwqSBe1Yq-FBUWVBOUdvaThEU1k
>>
>>
> I haven't my solution yet but I'm closer now, your suggestions were really helpful.
Marta

        [[alternative HTML version deleted]]

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

Re: Error in bbox(newdata) : object not a >= 2-column array

Tue, 03/14/2017 - 03:00
Looks to me as if you made a typo in your script, possibly with a quote,
but I'd need to see the full script, and preferably have access the
necessary data in order to run it, to be completely sure.

On 14/03/17 09:54, [hidden email] wrote:
>> Dear David,
> Dear R users,
> Now R calculates var1.pred as a constatnt:
> kriged=krige(Nmin~1,locationsDD,grid,model=fit)
>
>> kriged["var1.pred"]
>
>
> 925  (339003.4, 2988160)  11.22117
> 926  (339013.4, 2988160)  11.22117
> 927  (339023.4, 2988160)  11.22117
> 928  (339033.4, 2988160)  11.22117
> 929  (339043.4, 2988160)  11.22117
> 930  (339053.4, 2988160)  11.22117
> 931  (339063.4, 2988160)  11.22117
> 932  (339073.4, 2988160)  11.22117
> 933  (339083.4, 2988160)  11.22117
> ........................................
> 7482 (339293.4, 2989180)  11.22117
> 7483 (339303.4, 2989180)  11.22117
> 7484 (339313.4, 2989180)  11.22117
> 7485 (339323.4, 2989180)  11.22117
> 7486 (339333.4, 2989180)  11.22117
> 7487 (339343.4, 2989180)  11.22117
> 7488 (339353.4, 2989180)  11.22117
>>
> spplot(kriged["var1.pred"])
>  Error: unexpected string constant in:
> "
> spplot(kriged[""
>
> And it cann't give ssplot
> Please if anyone can help me!
>
>
> I was getting this error as well.  The cause was either not having the
>> coordinates and coordinate reference system (CRS) set or having the wrong
>> CRS.  I did both at the same time so I am not sure which it was.
>>
>> Is your 'grid' a SpatialPixels object?
>>
>> Did you run a line like 'coordinates(grid) <- ~x + y'  ?
>>
>> If so, did you then set the CRS?
>>
>>
>>
>> On Sat, Mar 11, 2017 at 7:38 AM, <[hidden email]> wrote:
>>
>>> Dear R-Users,
>>> Please anyone to help me with this error:
>>>
>>> kriged=krige(Nmin~1,locationsDD,grid,model=vgm(1,"Exp",150,1))
>>> Error in bbox(newdata) : object not a >= 2-column array
>>>
>>> Thank you so much in advance!
>>> Kind Regards
>>> Nina Philipova
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> [hidden email]
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>>
>>
>> --
>> David Warner
>> Research Fisheries Biologist
>> U.S.G.S. Great Lakes Science Center
>> 1451 Green Road
>> Ann Arbor, MI 48105
>> 734-214-9392
>>
>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo--
Edzer Pebesma
Institute for Geoinformatics  (ifgi),  University of Münster
Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
Journal of Statistical Software:   http://www.jstatsoft.org/
Computers & Geosciences:   http://elsevier.com/locate/cageo/


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
signature.asc (484 bytes) Download Attachment

Re: Error in bbox(newdata) : object not a >= 2-column array

Tue, 03/14/2017 - 02:54
> Dear David,
Dear R users,
Now R calculates var1.pred as a constatnt:
kriged=krige(Nmin~1,locationsDD,grid,model=fit)

> kriged["var1.pred"]


925  (339003.4, 2988160)  11.22117
926  (339013.4, 2988160)  11.22117
927  (339023.4, 2988160)  11.22117
928  (339033.4, 2988160)  11.22117
929  (339043.4, 2988160)  11.22117
930  (339053.4, 2988160)  11.22117
931  (339063.4, 2988160)  11.22117
932  (339073.4, 2988160)  11.22117
933  (339083.4, 2988160)  11.22117
........................................
7482 (339293.4, 2989180)  11.22117
7483 (339303.4, 2989180)  11.22117
7484 (339313.4, 2989180)  11.22117
7485 (339323.4, 2989180)  11.22117
7486 (339333.4, 2989180)  11.22117
7487 (339343.4, 2989180)  11.22117
7488 (339353.4, 2989180)  11.22117
>
spplot(kriged["var1.pred"])
 Error: unexpected string constant in:
"
spplot(kriged[""

And it cann't give ssplot
Please if anyone can help me!


I was getting this error as well.  The cause was either not having the
> coordinates and coordinate reference system (CRS) set or having the wrong
> CRS.  I did both at the same time so I am not sure which it was.
>
> Is your 'grid' a SpatialPixels object?
>
> Did you run a line like 'coordinates(grid) <- ~x + y'  ?
>
> If so, did you then set the CRS?
>
>
>
> On Sat, Mar 11, 2017 at 7:38 AM, <[hidden email]> wrote:
>
>> Dear R-Users,
>> Please anyone to help me with this error:
>>
>> kriged=krige(Nmin~1,locationsDD,grid,model=vgm(1,"Exp",150,1))
>> Error in bbox(newdata) : object not a >= 2-column array
>>
>> Thank you so much in advance!
>> Kind Regards
>> Nina Philipova
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>
>
> --
> David Warner
> Research Fisheries Biologist
> U.S.G.S. Great Lakes Science Center
> 1451 Green Road
> Ann Arbor, MI 48105
> 734-214-9392
>
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Variogram & Kriging.txt (3K) Download Attachment

Re: Avoid axis labels, ticks and marks in ggRGB() plot

Mon, 03/13/2017 - 07:25
Yes! Thanks!
I knew there had to be a better way. For my set of plots I can set
t <- theme(
  axis.text = element_blank(),
  axis.ticks = element_blank(),
  axis.title = element_blank())

and then

ggRGB(rlogo, r=1, g=2, b=3) + t +ggtitle("patatin")

Agus



On Mon, Mar 13, 2017 at 2:05 PM, Benjamin Leutner
<[hidden email]> wrote:
> The problem with your first approach is that there is no "data layer" which
> could span the coordinate axes.
> Since ggRGB by default uses an annotation_raster() in order to keep your
> "fill" dimension free for other data,  the coordinate limits are empty,
> hence the plotted raster is lost in the dimensionless void and you get see
> nothing ;-)
>
> I would recommend to start your plot with ggRGB like this:
>
> ggRGB(rlogo, r=1, g=2, b=3) +
>   theme(
>     axis.text = element_blank(),
>     axis.ticks = element_blank(),
>     axis.title = element_blank())
>
>
> If you do want to use it with ggLayer = TRUE, there should be at least one
> "geom" with actual data in it (or you'd have to set ggRGB(...,
> geom_raster=TRUE), but I don't see any benefit in that and it would give you
> a huge legend with one entry per rgb color...)
>
> best,
> Benjamin
>
>
> On 03/13/2017 12:10 PM, Agustin Lobo wrote:
>>
>> In order to save space, I prefer to avoid coordinates in the plot
>> generated by ggRGB. Is there any way? I've tried:
>>
>> ggplot() +
>> ggRGB(rlogo, r=1, g=2, b=3,ggLayer=TRUE) +
>>    theme(
>>      axis.text = element_blank(),
>>      axis.ticks = element_blank(),
>>      axis.title = element_blank())
>>
>> but get an empty plot.
>> I can circumvent the problem as follows, but it is kind of ugly and
>> complicated,
>> there must be another way:
>>
>> a <- as(extent(rlogo), 'SpatialPolygons')
>> y <- fortify(a,region=id)
>> ggplot(data=y,aes(y=lat, x=long)) +
>>    ggRGB(rlogo, r=1, g=2, b=3,ggLayer=TRUE) +
>>    geom_polygon(colour="yellow", fill=NA) +
>>    guides(fill=FALSE) +
>>    coord_fixed(ratio=1) +
>>    theme(
>>      axis.text = element_blank(),
>>      axis.ticks = element_blank(),
>>      axis.title = element_blank())
>>
>> Agus
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> --
> Benjamin Leutner M.Sc.
> PhD candidate
>
> Department of Remote Sensing
> University of Wuerzburg
> Campus Hubland Nord 86
> 97074 Wuerzburg, Germany
>
> Tel: +49-(0)931-31 89594
> Fax: +49-(0)931-31 89594-0
> Email: [hidden email]
> Web: http://www.fernerkundung.uni-wuerzburg.de
>
> _______________________________________________
> 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: Avoid axis labels, ticks and marks in ggRGB() plot

Mon, 03/13/2017 - 07:05
The problem with your first approach is that there is no "data layer"
which could span the coordinate axes.
Since ggRGB by default uses an annotation_raster() in order to keep your
"fill" dimension free for other data,  the coordinate limits are empty,
hence the plotted raster is lost in the dimensionless void and you get
see nothing ;-)

I would recommend to start your plot with ggRGB like this:

ggRGB(rlogo, r=1, g=2, b=3) +
   theme(
     axis.text = element_blank(),
     axis.ticks = element_blank(),
     axis.title = element_blank())


If you do want to use it with ggLayer = TRUE, there should be at least
one "geom" with actual data in it (or you'd have to set ggRGB(...,
geom_raster=TRUE), but I don't see any benefit in that and it would give
you a huge legend with one entry per rgb color...)

best,
Benjamin

On 03/13/2017 12:10 PM, Agustin Lobo wrote:
> In order to save space, I prefer to avoid coordinates in the plot
> generated by ggRGB. Is there any way? I've tried:
>
> ggplot() +
> ggRGB(rlogo, r=1, g=2, b=3,ggLayer=TRUE) +
>    theme(
>      axis.text = element_blank(),
>      axis.ticks = element_blank(),
>      axis.title = element_blank())
>
> but get an empty plot.
> I can circumvent the problem as follows, but it is kind of ugly and complicated,
> there must be another way:
>
> a <- as(extent(rlogo), 'SpatialPolygons')
> y <- fortify(a,region=id)
> ggplot(data=y,aes(y=lat, x=long)) +
>    ggRGB(rlogo, r=1, g=2, b=3,ggLayer=TRUE) +
>    geom_polygon(colour="yellow", fill=NA) +
>    guides(fill=FALSE) +
>    coord_fixed(ratio=1) +
>    theme(
>      axis.text = element_blank(),
>      axis.ticks = element_blank(),
>      axis.title = element_blank())
>
> Agus
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Benjamin Leutner M.Sc.
PhD candidate

Department of Remote Sensing
University of Wuerzburg
Campus Hubland Nord 86
97074 Wuerzburg, Germany

Tel: +49-(0)931-31 89594
Fax: +49-(0)931-31 89594-0
Email: [hidden email]
Web: http://www.fernerkundung.uni-wuerzburg.de

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

Avoid axis labels, ticks and marks in ggRGB() plot

Mon, 03/13/2017 - 05:10
In order to save space, I prefer to avoid coordinates in the plot
generated by ggRGB. Is there any way? I've tried:

ggplot() +
ggRGB(rlogo, r=1, g=2, b=3,ggLayer=TRUE) +
  theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    axis.title = element_blank())

but get an empty plot.
I can circumvent the problem as follows, but it is kind of ugly and complicated,
there must be another way:

a <- as(extent(rlogo), 'SpatialPolygons')
y <- fortify(a,region=id)
ggplot(data=y,aes(y=lat, x=long)) +
  ggRGB(rlogo, r=1, g=2, b=3,ggLayer=TRUE) +
  geom_polygon(colour="yellow", fill=NA) +
  guides(fill=FALSE) +
  coord_fixed(ratio=1) +
  theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    axis.title = element_blank())

Agus

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

Does ggmap work with sf objects through geom_sf?

Sun, 03/12/2017 - 07:29
Dear list members,

Does ggmap work with sf objects through geom_sf?

I got the following error:

nc <- st_read(system.file("shape/nc.shp", package="sf"))

nclocation <- c(-80, 34)

ncmap <- get_map(location = nclocation, zoom = 6)

ggmap(ncmap) + geom_sf(data = nc)

Error in eval(expr, envir, enclos) : object 'lon' not found


--
*Manuel Spínola, Ph.D.*
Instituto Internacional en Conservación y Manejo de Vida Silvestre
Universidad Nacional
Apartado 1350-3000
Heredia
COSTA RICA
[hidden email] <[hidden email]>
[hidden email]
Teléfono: (506) 8706 - 4662
Personal website: Lobito de río <https://sites.google.com/site/lobitoderio/>
Institutional website: ICOMVIS <http://www.icomvis.una.ac.cr/>

        [[alternative HTML version deleted]]

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

Re: Error in bbox(newdata) : object not a >= 2-column array

Sat, 03/11/2017 - 09:13
I was getting this error as well.  The cause was either not having the
coordinates and coordinate reference system (CRS) set or having the wrong
CRS.  I did both at the same time so I am not sure which it was.

Is your 'grid' a SpatialPixels object?

Did you run a line like 'coordinates(grid) <- ~x + y'  ?

If so, did you then set the CRS?



On Sat, Mar 11, 2017 at 7:38 AM, <[hidden email]> wrote:

> Dear R-Users,
> Please anyone to help me with this error:
>
> kriged=krige(Nmin~1,locationsDD,grid,model=vgm(1,"Exp",150,1))
> Error in bbox(newdata) : object not a >= 2-column array
>
> Thank you so much in advance!
> Kind Regards
> Nina Philipova
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>


--
David Warner
Research Fisheries Biologist
U.S.G.S. Great Lakes Science Center
1451 Green Road
Ann Arbor, MI 48105
734-214-9392

        [[alternative HTML version deleted]]

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

Re: Error in bbox(newdata) : object not a >= 2-column array

Sat, 03/11/2017 - 08:03
On Sat, Mar 11, 2017 at 12:38 PM,  <[hidden email]> wrote:
> Dear R-Users,
> Please anyone to help me with this error:
>
> kriged=krige(Nmin~1,locationsDD,grid,model=vgm(1,"Exp",150,1))
> Error in bbox(newdata) : object not a >= 2-column array
>
> Thank you so much in advance!

 Hard to debug without your data. The error message is telling us that
something is wrong with your input, but we can't see your data....

 If you can't share your data, at least share what you get from
summary(grid) and summary(locationsDD).

 You should also confirm what package `krige` is coming from - I
assume its gstat.

If you can get back to the mailing list with this information, we
might be able to help, otherwise we're guessing!

Barry

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

Error in bbox(newdata) : object not a >= 2-column array

Sat, 03/11/2017 - 06:38
Dear R-Users,
Please anyone to help me with this error:

kriged=krige(Nmin~1,locationsDD,grid,model=vgm(1,"Exp",150,1))
Error in bbox(newdata) : object not a >= 2-column array

Thank you so much in advance!
Kind Regards
Nina Philipova
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
variogram2.txt (2K) Download Attachment

Re: Global and Local Moran: how can I calculate them with different spatial threshold?

Fri, 03/10/2017 - 17:45
Perfect, thank you very much to everybody,
Next week I will work on my data. Now everything is much clearer!

Regards

Il giorno ven 10 mar 2017 alle 21:44 Roger Bivand <[hidden email]> ha
scritto:

> On Fri, 10 Mar 2017, Maurizio Marchi wrote:
>
> > Hello everybody,
> > I need to calculate the Global and Local Moran indices in R using a
> > variable distance as threshold: I have a SpatialPointsDataFrame with
> > almost 300 points and I want to calculate the Global Moran index using
> > 4 different distances (e.g 5 km - 10 km - 50 km - 100 km).
> > I know the 'ape' and 'spdep' packages but it seems that no adjustment
> > can be done concerning the spatial width to be considered...
>
> Please examine ?spdep::dnearneigh, especially d1= and d2=, which do
> exactly what you ask for. The functions in spdep are modularised, first
> construct the neighbour object, then the weights list object, then the
> Moran tests. Note that Moran tests should really be run on regression
> residuals (lm.morantest(), localmoran.sad() or localmoran.exact()), and
> remember to adjust p.values for banded local tests (many tests using the
> same data affect tabulated significance levels). The ape function makes
> undocumented assumptions about what you may want, spdep requires that you
> know what you want to do.
>
> Hope this clarifies,
>
> Roger
>
> > Thanks
> >
> >
>
> --
> 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]
> Editor-in-Chief of The R Journal, https://journal.r-project.org/index.html
> http://orcid.org/0000-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
> http://depsy.org/person/444584
> --
Maurizio Marchi
Skype ID: maurizioxyz
linux user 552742

        [[alternative HTML version deleted]]

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

Re: Global and Local Moran: how can I calculate them with different spatial threshold?

Fri, 03/10/2017 - 14:44
On Fri, 10 Mar 2017, Maurizio Marchi wrote:

> Hello everybody,
> I need to calculate the Global and Local Moran indices in R using a
> variable distance as threshold: I have a SpatialPointsDataFrame with
> almost 300 points and I want to calculate the Global Moran index using
> 4 different distances (e.g 5 km - 10 km - 50 km - 100 km).
> I know the 'ape' and 'spdep' packages but it seems that no adjustment
> can be done concerning the spatial width to be considered...

Please examine ?spdep::dnearneigh, especially d1= and d2=, which do
exactly what you ask for. The functions in spdep are modularised, first
construct the neighbour object, then the weights list object, then the
Moran tests. Note that Moran tests should really be run on regression
residuals (lm.morantest(), localmoran.sad() or localmoran.exact()), and
remember to adjust p.values for banded local tests (many tests using the
same data affect tabulated significance levels). The ape function makes
undocumented assumptions about what you may want, spdep requires that you
know what you want to do.

Hope this clarifies,

Roger

> Thanks
>
>

--
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]
Editor-in-Chief of The R Journal, https://journal.r-project.org/index.html
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
http://depsy.org/person/444584

_______________________________________________
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

Global and Local Moran: how can I calculate them with different spatial threshold?

Fri, 03/10/2017 - 12:16
Hello everybody,
I need to calculate the Global and Local Moran indices in R using a
variable distance as threshold: I have a SpatialPointsDataFrame with
almost 300 points and I want to calculate the Global Moran index using
4 different distances (e.g 5 km - 10 km - 50 km - 100 km).
I know the 'ape' and 'spdep' packages but it seems that no adjustment
can be done concerning the spatial width to be considered...
Thanks

--
Maurizio Marchi
Skype ID: maurizioxyz
linux user 552742

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

Re: noob question: access single region from .shp file

Thu, 03/09/2017 - 21:22
On Thu, Mar 9, 2017 at 6:00 AM, <[hidden email]> wrote:
>
> Message: 2
> Date: Wed, 08 Mar 2017 15:45:45 +0200
> From: Brandon Payne <[hidden email]>
> To: [hidden email]
> Subject: [R-sig-Geo] noob question: access single region from .shp
>         file
> Message-ID: <[hidden email]>
> Content-Type: text/plain
>
> How do you access a single region in a shape file?
>
See changes below.

Kent

## libraries
> library(maptools)
> library(RColorBrewer)
> colors <- brewer.pal(9, "BuGn")
> library(ggmap)
>
> ## google maps
> mapImageData3 <- get_map(location = c(lon = 35.1660235,
>                                       lat = 31.32226),
>                          color="color",
>                          source = "google",
>                          maptype = "satellite",
>                          zoom = 7)
>
> ## read shape file
> IsraelPolygon <- readShapePoly("./includes/ISR_adm/ISR_adm1.shp")
> IsraelPoints  <- fortify(IsraelPolygon)
>
IsraelPoints <- fortify(IsraelPolygon, , region='NAME_1')

## plot shapes on top of map, output to graphics device
> ggmap(mapImageData3,
>     extent = "device",
>     ylab = "Latitude",
>     xlab = "Longitude") +
>   geom_polygon(aes(x=long,
>                    y=lat,
>                    group=group),
>                data=IsraelPoints,
>                color=colors[9],
>                fil=  colors[6],
>
fill=  ifelse(IsraelPoints$id=='Golan', 'red', colors[6]),


>                alpha=0.5) +
>   labs(x="Longitude",
>        y="Latitude"

        [[alternative HTML version deleted]]

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

Pages