Re: [R-sig-Geo] [FORGED] different models

2016-04-06 Thread Rolf Turner 

On 06/04/16 22:00, Virginia Morera Pujol wrote:

Hi all,

This might be a very dumb question that shows I have very little idea of
what I am talking about, but I'll risk it:

What is the difference between fitting a model using these 3 different
syntaxes?

1/ fit1 <- ppm(ppp, ~covariate),
2/ fit2 <- ppm(ppp, ~x+y+Z, covariates=list(Z=covariate))
3/ fit3 <- ppm(ppp, ~x+y+covariate)

where ppp is my point pattern and "covariate" is a pixel image? I realise
the outputs of 2 and 3 are the same and different to that of 1, so I guess
the question really is

a/ Is there any difference, practical or in the actual computations of the
model, between using 2 and 3?
b/ What is the difference between (2&3) and 1?


(1) There is essentially no difference between fits 2 & 3.  The fit 2 
syntax is provided so that the user can have the relevant covariates 
bundled up in a list without any need to extract these covariates from 
that list.   With the fit 2 syntax you don't need to have all covariates 
present in your workspace.


E.g.: fit <- ppm(bei ~ elev + grad, data=bei.extra)

(2) The fit 2 syntax is essentially the same as that used by lm() and 
glm() and was designed in imitation thereof.


(3) The preferred structure of a call to ppm() is

fit2 <- ppm(ppp ~ x + y + Z, data=list(Z=covariate))

Note:  "data" rather than "covariates"; no comma between the name of the 
response variable ("ppp") and the formula.


This makes the syntax identical to that of lm() and glm().

The syntax that you used is a remnant of earlier versions of spatstat 
and remains acceptable for reasons of backward compatibility.


(4) The difference between model 1 and models 2 and 3 is that models 2 
and 3 involve the Cartesian coordinates "x" and "y".  Model 1 is such 
that the model intensity takes the form


   exp(beta_0 + beta_1 * covariate)

In models 2 and 3 the model intensity takes the (more complex) form

   exp(beta_0 + beta_1 * x + beta_2 *y beta_3 * covariate)

Note that "x" and "y" are *reserved* names.  You cannot use these names 
for any covariates *other than* the Cartesian coordinates.


(5) The name "covariate" is probably *not* a good name for a covariate.
As fortune(77) puts it "Would you call your dog 'dog'?"

(6) Likewise (and even more so) "ppp" is *not* a good name for a point 
pattern, since it clashes the name of the creator function ppp().


cheers,

Rolf Turner

--
Technical Editor ANZJS
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] readOGR not reading in prj file?

2016-04-06 Thread Vinh Nguyen 
On Wed, Apr 6, 2016 at 10:42 AM, Edzer Pebesma
 wrote:
> When I replace the .prj file of an arbitrary shapefile with these
> contents, it gets read well. Have you tried with current R (3.2.4) and
> rgdal (1.1.7)?


Thanks Edzer, upgrading rgdal to 1.1.8 works.  Consider this issue close.

-- Vinh

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] readOGR not reading in prj file?

2016-04-06 Thread Edzer Pebesma 


On 06/04/16 19:33, Vinh Nguyen wrote:
> 
> On Apr 6, 2016 10:27 AM, "Edzer Pebesma"  > wrote:
>>
>> What are the contents of the 2013Spring_Flood_Drawn.prj file?
>>
> 
> Here's the content of the prj file:
> PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]
> 
> Thanks so much for your help.
> 

When I replace the .prj file of an arbitrary shapefile with these
contents, it gets read well. Have you tried with current R (3.2.4) and
rgdal (1.1.7)?
-- 
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/
Spatial Statistics Society http://www.spatialstatistics.info



signature.asc
Description: OpenPGP digital signature
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] readOGR not reading in prj file?

2016-04-06 Thread Vinh Nguyen 
On Apr 6, 2016 10:27 AM, "Edzer Pebesma" 
wrote:
>
> What are the contents of the 2013Spring_Flood_Drawn.prj file?
>

Here's the content of the prj file:
PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]

Thanks so much for your help.

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] readOGR not reading in prj file?

2016-04-06 Thread Vinh Nguyen 
On Wed, Apr 6, 2016 at 8:58 AM, Denis mutiibwa  wrote:
> readOGR(dir.path, "argv2")
> it seems "argv2" has to be the shapefile name not a variable name


My 'floodshapedir' and 'shapebase' are both variables with path stored
in.  Here's the function call with the strings explicitly passed.  I
don't think it's an issue with these paths.

> dFloodRisk <- readOGR('FloodShapes', '2013Spring_Flood_Drawn', verbose=TRUE)
OGR data source with driver: ESRI Shapefile
Source: "FloodShapes", layer: "2013Spring_Flood_Drawn"
with 8 features
It has 2 fields
> proj4string(dFloodRisk)
[1] NA

-- Vinh

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] readOGR not reading in prj file?

2016-04-06 Thread Denis mutiibwa via R-sig-Geo 
readOGR(dir.path, "argv2")
it seems "argv2" has to be the shapefile name not a variable name   Denis 
Mutiibwa 

On Wednesday, April 6, 2016 11:54 AM, Vinh Nguyen  
wrote:
 

 On Wed, Apr 6, 2016 at 8:41 AM, Denis mutiibwa  wrote:
>> Where is "shapebase" coming from?
>> If those are the files in the directory, your scripts should as below:
>> dir.path <- directory path
>> dFloodRisk <- readOGR(dir.path, "2014Spring_Flood_Drawn")
>> proj4string(dFloodRisk)
>
>I'm reading the files in exactly as you are suggesting (floodshapedir
>and shapebase are variables) as seen in the output below:
>
>> dFloodRisk <- readOGR(floodshapedir, shapebase)
>OGR data source with driver: ESRI Shapefile
>Source: "FloodShapes", layer: "2013Spring_Flood_Drawn"
>
>Just not sure why the prj file is not being read in.  Thanks.>
>
>-- Vinh


  
[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] readOGR not reading in prj file?

2016-04-06 Thread Vinh Nguyen 
On Wed, Apr 6, 2016 at 8:41 AM, Denis mutiibwa  wrote:
> Where is "shapebase" coming from?
> If those are the files in the directory, your scripts should as below:
> dir.path <- directory path
> dFloodRisk <- readOGR(dir.path, "2014Spring_Flood_Drawn")
> proj4string(dFloodRisk)

I'm reading the files in exactly as you are suggesting (floodshapedir
and shapebase are variables) as seen in the output below:

> dFloodRisk <- readOGR(floodshapedir, shapebase)
OGR data source with driver: ESRI Shapefile
Source: "FloodShapes", layer: "2013Spring_Flood_Drawn"

Just not sure why the prj file is not being read in.  Thanks.

-- Vinh

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Faster way to extract raster values using multiple polygons?

2016-04-06 Thread Ben Tupper 
Hi,

The resolution of your raster data (1 degree) is much more coarse than what 
your polygons represent.  Could you short-circuit the process by assuming that 
the temp at the centroid of each polygon would suitably represent the mean 
temperature across each polygon?  Unless you have some much bigger polygons, I 
can't imagine it will be very far off. If so, then you could pretty quickly 
extract the values for each layer in the raster at each centroid.  Perhaps like 
this?

cents <- coordinates(br_sub)
v <- extract(b, cents)

Is that close enough?

Cheers,
Ben

> On Apr 6, 2016, at 8:00 AM, Thiago V. dos Santos via R-sig-Geo 
>  wrote:
> 
> Dear all,
> 
> I am trying to extract a time series from a raster brick using multiple 
> polygons.
> 
> The raster brick is a global temperature netcdf file with 1200 layers 
> (months) and the polygons are each of the 5567 municipalities in Brazil. I 
> intend to extract the temperature time series for each municipality.
> 
> As a final result, I would like to have a data frame containing: -the name 
> and coordinates of the municipality, -the date tag from the raster, and -the 
> extracted temperature value.
> 
> My initial, VERY SLOW attempt was something like this:
> 
> 
> library(raster)
> 
> # Create some sample raster data
> idx <- seq(as.Date("1960/1/1"), as.Date("1990/12/01"), by = "month")
> r <- raster(ncol=360, nrow=180)
> b <- brick(lapply(1:length(idx), function(x) setValues(r, runif(ncell(r)
> b <- setZ(b, idx)
> 
> # Import polygon data
> br <- getData("GADM", country="Brazil", level=2) # about 7MB to download
> 
> # Subset the SMALLEST state - only 75 municipalities
> br_sub <- br[br$NAME_1=="Sergipe" & !is.na(br$NAME_1),]
> plot(br_sub)
> 
> # Now let's extract the data for each municipality
> beginCluster()
> e <- extract(b, br_sub, fun=mean, df=T)
> endCluster()
> 
> 
> Even using the smallest state possible, this example takes about 20 minutes 
> to run on a dual-core 2.5GHz Macbook Pro using four threads. As a reference, 
> there are states with over 850 municipalities. And remember, in total there 
> are over 5500 municipalities I need to extract the data for.
> 
> I have the feeling that my code is not very optimized.
> 
> Has anybody ever dealt with this amount of data in this kind of raster 
> operation? Is there any fastest way to achieve my expected result?
> Thanks in advance,
> -- Thiago V. dos Santos
> 
> PhD student
> Land and Atmospheric Science
> 
> University of Minnesota
> 
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> 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
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] readOGR not reading in prj file?

2016-04-06 Thread Denis mutiibwa via R-sig-Geo 
Where is "shapebase" coming from?If those are the files in the directory, your 
scripts should as below:dir.path <- directory pathdFloodRisk <- 
readOGR(dir.path, "2014Spring_Flood_Drawn")proj4string(dFloodRisk) 
Denis Mutiibwa 

On Wednesday, April 6, 2016 11:32 AM, Vinh Nguyen  
wrote:
 

 On Wed, Apr 6, 2016 at 8:24 AM, Denis mutiibwa  wrote:
>> Hi Vinh,
>> proj4string function is giving you NA because "shapebase" shapefile doesn't
>> have an auxiliary file "shapebase.prj"
>>
>> Denis Mutiibwa,
>
>
>Thanks for your response Denis.  The following files are all in the directory:
>2014Spring_Flood_Drawn.cpg  2014Spring_Flood_Drawn.sbx
>2014Spring_Flood_Drawn.dbf  2014Spring_Flood_Drawn.shp
>2014Spring_Flood_Drawn.prj  2014Spring_Flood_Drawn.shp.xml
>2014Spring_Flood_Drawn.sbn  2014Spring_Flood_Drawn.shx
>
>So, the prj file is there...or am I misunderstanding?  Thanks.>
>
>-- Vinh


  
[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] readOGR not reading in prj file?

2016-04-06 Thread Vinh Nguyen 
On Wed, Apr 6, 2016 at 8:24 AM, Denis mutiibwa  wrote:
> Hi Vinh,
> proj4string function is giving you NA because "shapebase" shapefile doesn't
> have an auxiliary file "shapebase.prj"
>
> Denis Mutiibwa,


Thanks for your response Denis.  The following files are all in the directory:
2014Spring_Flood_Drawn.cpg  2014Spring_Flood_Drawn.sbx
2014Spring_Flood_Drawn.dbf  2014Spring_Flood_Drawn.shp
2014Spring_Flood_Drawn.prj  2014Spring_Flood_Drawn.shp.xml
2014Spring_Flood_Drawn.sbn  2014Spring_Flood_Drawn.shx

So, the prj file is there...or am I misunderstanding?  Thanks.

-- Vinh

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] readOGR not reading in prj file?

2016-04-06 Thread Denis mutiibwa via R-sig-Geo 
Hi Vinh,proj4string function is giving you NA because "shapebase" shapefile 
doesn't have an auxiliary file "shapebase.prj" Denis Mutiibwa, 

On Wednesday, April 6, 2016 10:46 AM, Vinh Nguyen  
wrote:
 

 >Hi,

>I'm currently using rgdal v 1.1-1 on R 3.2.2.  It appears my
>projection file is not being read in by readOGR:
>
>> dFloodRisk <- readOGR(floodshapedir, shapebase)
>OGR data source with driver: ESRI Shapefile
>Source: "FloodShapes", layer: "2013Spring_Flood_Drawn"
>with 8 features
>It has 2 fields
> proj4string(dFloodRisk)
>[1] NA
>
>I have to manually pass in the CRS:
>> dFloodRisk <- readOGR(floodshapedir, paste0(curYear, 'Spring_Flood_Drawn'), 
>> p4s='+init=EPSG:3857')
>OGR data source with driver: ESRI Shapefile
>Source: "FloodShapes", layer: "2013Spring_Flood_Drawn"
>with 8 features
>It has 2 fields
>> proj4string(dFloodRisk)
>[1] "+init=EPSG:3857 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0
>+lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"
>
>Any thoughts on why readOGR doesn't automatically read in the prj
>file?  Here's the content of the prj file:
>PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Gree>nwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER[">>Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]
>
>Thanks so much for your help.
>
>-- Vinh

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo


  
[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] ERA Interim Daily Data - Query historical data for a point (defined by lat and lon)

2016-04-06 Thread Sébastien Piguet 
Hi,

I am currently trying to query ECMWF data (
http://apps.ecmwf.int/datasets/data/interim-full-daily/levtype=sfc/) for a
point (lon,lat) ~ (-71.5333,-32.7723). I query an area around this point as
showed below:
[image: Inline image 1]
Then, I use R to open the grip file:

library(rgdal)
library(maptools)
library(raster)
library(RCurl)
library(stringr)

setwd("XX")
#Initialize raster stack
rs<-stack()

#Get filenames in the folder
filePath<-"YYY"
fileNms<-list.files(filePath)

#Define Point of interest lon / lat
point <- cbind(-71.5333,-32.7723)

# DB loading
grib <- readGDAL(paste(filePath,fileNms[1], sep=""))

# Selection of significant height of total swell = param 237 (
http://www.ecmwf.int/sites/default/files/elibrary/2011/8174-era-interim-archive-version-20.pdf
)
rLayer<-raster(grib, 237)
rLayer


I get the following results:

> summary(rLayer)  band62
Min.10.25290
1st Qu. 10.33603
Median  10.41916
3rd Qu. 10.41916
Max.10.41916
NA's15.0> summary(rasterToPoints(rLayer))   x
 y  band62
 Min.   :288.0   Min.   :-33.75   Min.   :10.25
 1st Qu.:288.0   1st Qu.:-33.00   1st Qu.:10.34
 Median :288.0   Median :-32.25   Median :10.42
 Mean   :288.2   Mean   :-32.75   Mean   :10.36
 3rd Qu.:288.4   3rd Qu.:-32.25   3rd Qu.:10.42
 Max.   :288.8   Max.   :-32.25   Max.   :10.42

My problems are the following:
1) Latitude is similar from what I expected  but not the longitude.

2) Data about date does not appear

Thank you so much for your help !

Sebastien
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] can't read in shapefile

2016-04-06 Thread Vinh Nguyen 
Hi,

I have two shapefiles, both generated from ArcGIS using similar
methods.  I was able to read in one file using readOGR but not the
other.  Any thoughts on how best to figure out why I can't read in the
file?  Here's the error:

> dFloodRisk <- readOGR(floodshapedir, shapefilebase, verbose=TRUE)
Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding,
use_iconv = use_iconv,  :
  Cannot open layer
> traceback()
5: .Call("ogrInfo", as.character(dsn), as.character(layer), PACKAGE = "rgdal")
4: ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv = use_iconv,
   swapAxisOrder = swapAxisOrder, require_geomType = require_geomType)
3: withCallingHandlers(expr, message = function(c)
invokeRestart("muffleMessage"))
2: suppressMessages(ogr_info <- ogrInfo(dsn = dsn, layer = layer,
   encoding = encoding, use_iconv = use_iconv, swapAxisOrder =
swapAxisOrder,
   require_geomType = require_geomType))
1: readOGR(floodshapedir, shapefilebase,
   verbose = TRUE)

Thanks for your help.

-- Vinh

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] readOGR not reading in prj file?

2016-04-06 Thread Vinh Nguyen 
Hi,

I'm currently using rgdal v 1.1-1 on R 3.2.2.  It appears my
projection file is not being read in by readOGR:

> dFloodRisk <- readOGR(floodshapedir, shapebase)
OGR data source with driver: ESRI Shapefile
Source: "FloodShapes", layer: "2013Spring_Flood_Drawn"
with 8 features
It has 2 fields
> proj4string(dFloodRisk)
[1] NA

I have to manually pass in the CRS:
> dFloodRisk <- readOGR(floodshapedir, paste0(curYear, 'Spring_Flood_Drawn'), 
> p4s='+init=EPSG:3857')
OGR data source with driver: ESRI Shapefile
Source: "FloodShapes", layer: "2013Spring_Flood_Drawn"
with 8 features
It has 2 fields
> proj4string(dFloodRisk)
[1] "+init=EPSG:3857 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0
+lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"

Any thoughts on why readOGR doesn't automatically read in the prj
file?  Here's the content of the prj file:
PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]

Thanks so much for your help.

-- Vinh

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] problem with rasterToPolygons x worldclim

2016-04-06 Thread Michael Sumner 
On Thu, 7 Apr 2016 at 00:33 Michael Sumner  wrote:

> On Thu, 7 Apr 2016 at 00:22 Michael Sumner  wrote:
>
>> On Wed, 6 Apr 2016 at 23:45 Karla Shikev  wrote:
>>
>>>
>>> Hi Michael,
>>>
>>> Thanks for the tips. I tried the help for rasterToPolygons, but none of
>>> the options (e.g. dissolve=TRUE) made any difference. I tried gPolygonize
>>> and it worked, except for (as you predicted), if the edges meet the raster
>>> extents - but that definitely was an advance!
>>>
>>>
>> I see the issue, the fun argument to rasterToPolygons is masking out
>> values from that interval, but it's not setting them all to the same value
>> - so you still get individual pixel polygons unless you mask on presence in
>> the interval or not (as a binary):
>>
>> library(raster)
>>
>> r<-getData('worldclim', var='bio', res=10)$bio1
>>
>> e <- extent(-140,-100, 50, 60)
>>
>> xx <- crop(r,e)
>> ## mask out pixel values first
>> xx <- xx > 40 & xx < 60
>>
>> ## 576 polygons
>> slice1<- rasterToPolygons(xx, fun = function(x) {x == 1})
>> ## 1 polygon
>> slice2 <- rasterToPolygons(xx,  fun = function(x) {x == 1}, dissolve =
>> TRUE)
>>
>>
>> Still it's not a very nice polygon, if I can I'll try a different way.
>>
>> Cheers, Mike.
>>
>>
>>
>>
>
> You can use tricks to fill the outer edge so you get a valid contour line
> all the way around (it's not always going to work - see here for the
> explanation by whuber:
> http://gis.stackexchange.com/questions/61550/colouring-areas-between-vector-contours
> )
>
>

Sorry, forgot the code!

library(raster)

r<-getData('worldclim', var='bio', res=10)$bio1

e <- extent(-140,-100, 50, 60)

xx <- crop(r,e)

## buffer out a little
xx2 <- extend(xx, extent(xx) + res(xx) * 4, value = NA_real_)

## set missing to less than interval
xx2[is.na(xx2)] <- cellStats(xx2, min)
## now contour
cl <- rasterToContour(xx2 > 40 & xx2 < 60, level = 1)

pp <- rgeos::gUnionCascaded(rgeos::gPolygonize(cl))
pp
plot(pp, col = "grey")




> Is this better? Maybe - there are lots of other ways - maybe a buffer on
> this gives roughly what you need. Maybe you only need a total area so
> summing pixels after filtering is better?
>
> Would be nice to have contouring that reliably produced polygons, but it's
> a tricky problem. Other kinds of shape-finding might be better, with
> alphahull on the pixel points perhaps, but you quickly get into "geometry
> fudger" territory here rather than objective measures.
>
> HTH
>
>
>
>
>>
>>>
>>> I really appreciate your help and if you just point the way I can try to
>>> go after what needs to be done, but right now I´m stumped.
>>>
>>> best, Karla.
>>>
>>> On Mon, Apr 4, 2016 at 5:38 PM, Michael Sumner 
>>> wrote:
>>>
 Read the help for rasterToPolygons, alternatively try chaining
 rasterToContour and rgeos::gPolygonize.

 The latter may need coaxing if your edges meet the raster extents.

 Cheers, Mike

 On Tue, 5 Apr 2016, 05:32 Karla Shikev  wrote:

> Dear all,
>
> Here's an issue that is related to the previous one. In the commands
> below
> I'm trying to make a polygon for all the regions within a range of mean
> annual temps. However, rasterToPolygons will draw a square around each
> value in the original raster, rather than providing me the actual
> polygon,
> given that resolution. Any hints? Again, any help will be greatly
> appreciated.
>
> Karla
>
> _
>
> library(raster)
>
> r<-getData('worldclim', var='bio', res=10)$bio1
>
> e <- extent(-140,-100, 50, 60)
>
> xx<-crop(r,e)
>
> slice1<-rasterToPolygons(xx, fun = function(x) {x > 40 & x < 60})
>
> plot(slice1)
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> 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


>>> --
>> Dr. Michael Sumner
>> Software and Database Engineer
>> Australian Antarctic Division
>> 203 Channel Highway
>> Kingston Tasmania 7050 Australia
>>
>> --
> Dr. Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> 203 Channel Highway
> Kingston Tasmania 7050 Australia
>
> --
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
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] problem with rasterToPolygons x worldclim

2016-04-06 Thread Michael Sumner 
On Thu, 7 Apr 2016 at 00:22 Michael Sumner  wrote:

> On Wed, 6 Apr 2016 at 23:45 Karla Shikev  wrote:
>
>>
>> Hi Michael,
>>
>> Thanks for the tips. I tried the help for rasterToPolygons, but none of
>> the options (e.g. dissolve=TRUE) made any difference. I tried gPolygonize
>> and it worked, except for (as you predicted), if the edges meet the raster
>> extents - but that definitely was an advance!
>>
>>
> I see the issue, the fun argument to rasterToPolygons is masking out
> values from that interval, but it's not setting them all to the same value
> - so you still get individual pixel polygons unless you mask on presence in
> the interval or not (as a binary):
>
> library(raster)
>
> r<-getData('worldclim', var='bio', res=10)$bio1
>
> e <- extent(-140,-100, 50, 60)
>
> xx <- crop(r,e)
> ## mask out pixel values first
> xx <- xx > 40 & xx < 60
>
> ## 576 polygons
> slice1<- rasterToPolygons(xx, fun = function(x) {x == 1})
> ## 1 polygon
> slice2 <- rasterToPolygons(xx,  fun = function(x) {x == 1}, dissolve =
> TRUE)
>
>
> Still it's not a very nice polygon, if I can I'll try a different way.
>
> Cheers, Mike.
>
>
>
>

You can use tricks to fill the outer edge so you get a valid contour line
all the way around (it's not always going to work - see here for the
explanation by whuber:
http://gis.stackexchange.com/questions/61550/colouring-areas-between-vector-contours
)

Is this better? Maybe - there are lots of other ways - maybe a buffer on
this gives roughly what you need. Maybe you only need a total area so
summing pixels after filtering is better?

Would be nice to have contouring that reliably produced polygons, but it's
a tricky problem. Other kinds of shape-finding might be better, with
alphahull on the pixel points perhaps, but you quickly get into "geometry
fudger" territory here rather than objective measures.

HTH




>
>>
>> I really appreciate your help and if you just point the way I can try to
>> go after what needs to be done, but right now I´m stumped.
>>
>> best, Karla.
>>
>> On Mon, Apr 4, 2016 at 5:38 PM, Michael Sumner 
>> wrote:
>>
>>> Read the help for rasterToPolygons, alternatively try chaining
>>> rasterToContour and rgeos::gPolygonize.
>>>
>>> The latter may need coaxing if your edges meet the raster extents.
>>>
>>> Cheers, Mike
>>>
>>> On Tue, 5 Apr 2016, 05:32 Karla Shikev  wrote:
>>>
 Dear all,

 Here's an issue that is related to the previous one. In the commands
 below
 I'm trying to make a polygon for all the regions within a range of mean
 annual temps. However, rasterToPolygons will draw a square around each
 value in the original raster, rather than providing me the actual
 polygon,
 given that resolution. Any hints? Again, any help will be greatly
 appreciated.

 Karla

 _

 library(raster)

 r<-getData('worldclim', var='bio', res=10)$bio1

 e <- extent(-140,-100, 50, 60)

 xx<-crop(r,e)

 slice1<-rasterToPolygons(xx, fun = function(x) {x > 40 & x < 60})

 plot(slice1)

 [[alternative HTML version deleted]]

 ___
 R-sig-Geo mailing list
 R-sig-Geo@r-project.org
 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
>>>
>>>
>> --
> Dr. Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> 203 Channel Highway
> Kingston Tasmania 7050 Australia
>
> --
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
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] problem with rasterToPolygons x worldclim

2016-04-06 Thread Michael Sumner 
On Wed, 6 Apr 2016 at 23:45 Karla Shikev  wrote:

>
> Hi Michael,
>
> Thanks for the tips. I tried the help for rasterToPolygons, but none of
> the options (e.g. dissolve=TRUE) made any difference. I tried gPolygonize
> and it worked, except for (as you predicted), if the edges meet the raster
> extents - but that definitely was an advance!
>
>
I see the issue, the fun argument to rasterToPolygons is masking out values
from that interval, but it's not setting them all to the same value - so
you still get individual pixel polygons unless you mask on presence in the
interval or not (as a binary):

library(raster)

r<-getData('worldclim', var='bio', res=10)$bio1

e <- extent(-140,-100, 50, 60)

xx <- crop(r,e)
## mask out pixel values first
xx <- xx > 40 & xx < 60

## 576 polygons
slice1<- rasterToPolygons(xx, fun = function(x) {x == 1})
## 1 polygon
slice2 <- rasterToPolygons(xx,  fun = function(x) {x == 1}, dissolve = TRUE)


Still it's not a very nice polygon, if I can I'll try a different way.

Cheers, Mike.




>
> I really appreciate your help and if you just point the way I can try to
> go after what needs to be done, but right now I´m stumped.
>
> best, Karla.
>
> On Mon, Apr 4, 2016 at 5:38 PM, Michael Sumner  wrote:
>
>> Read the help for rasterToPolygons, alternatively try chaining
>> rasterToContour and rgeos::gPolygonize.
>>
>> The latter may need coaxing if your edges meet the raster extents.
>>
>> Cheers, Mike
>>
>> On Tue, 5 Apr 2016, 05:32 Karla Shikev  wrote:
>>
>>> Dear all,
>>>
>>> Here's an issue that is related to the previous one. In the commands
>>> below
>>> I'm trying to make a polygon for all the regions within a range of mean
>>> annual temps. However, rasterToPolygons will draw a square around each
>>> value in the original raster, rather than providing me the actual
>>> polygon,
>>> given that resolution. Any hints? Again, any help will be greatly
>>> appreciated.
>>>
>>> Karla
>>>
>>> _
>>>
>>> library(raster)
>>>
>>> r<-getData('worldclim', var='bio', res=10)$bio1
>>>
>>> e <- extent(-140,-100, 50, 60)
>>>
>>> xx<-crop(r,e)
>>>
>>> slice1<-rasterToPolygons(xx, fun = function(x) {x > 40 & x < 60})
>>>
>>> plot(slice1)
>>>
>>> [[alternative HTML version deleted]]
>>>
>>> ___
>>> R-sig-Geo mailing list
>>> R-sig-Geo@r-project.org
>>> 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
>>
>>
> --
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
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] problem with rasterToPolygons x worldclim

2016-04-06 Thread Karla Shikev 
Hi Michael,

Thanks for the tips. I tried the help for rasterToPolygons, but none of the
options (e.g. dissolve=TRUE) made any difference. I tried gPolygonize and
it worked, except for (as you predicted), if the edges meet the raster
extents - but that definitely was an advance!

I really appreciate your help and if you just point the way I can try to go
after what needs to be done, but right now I´m stumped.

best, Karla.

On Mon, Apr 4, 2016 at 5:38 PM, Michael Sumner  wrote:

> Read the help for rasterToPolygons, alternatively try chaining
> rasterToContour and rgeos::gPolygonize.
>
> The latter may need coaxing if your edges meet the raster extents.
>
> Cheers, Mike
>
> On Tue, 5 Apr 2016, 05:32 Karla Shikev  wrote:
>
>> Dear all,
>>
>> Here's an issue that is related to the previous one. In the commands below
>> I'm trying to make a polygon for all the regions within a range of mean
>> annual temps. However, rasterToPolygons will draw a square around each
>> value in the original raster, rather than providing me the actual polygon,
>> given that resolution. Any hints? Again, any help will be greatly
>> appreciated.
>>
>> Karla
>>
>> _
>>
>> library(raster)
>>
>> r<-getData('worldclim', var='bio', res=10)$bio1
>>
>> e <- extent(-140,-100, 50, 60)
>>
>> xx<-crop(r,e)
>>
>> slice1<-rasterToPolygons(xx, fun = function(x) {x > 40 & x < 60})
>>
>> plot(slice1)
>>
>> [[alternative HTML version deleted]]
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> 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
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Faster way to extract raster values using multiple polygons?

2016-04-06 Thread Thiago V. dos Santos via R-sig-Geo 
Dear all,

I am trying to extract a time series from a raster brick using multiple 
polygons.

The raster brick is a global temperature netcdf file with 1200 layers (months) 
and the polygons are each of the 5567 municipalities in Brazil. I intend to 
extract the temperature time series for each municipality.

As a final result, I would like to have a data frame containing: -the name and 
coordinates of the municipality, -the date tag from the raster, and -the 
extracted temperature value.

My initial, VERY SLOW attempt was something like this:


library(raster)

# Create some sample raster data
idx <- seq(as.Date("1960/1/1"), as.Date("1990/12/01"), by = "month")
r <- raster(ncol=360, nrow=180)
b <- brick(lapply(1:length(idx), function(x) setValues(r, runif(ncell(r)
b <- setZ(b, idx)

# Import polygon data
br <- getData("GADM", country="Brazil", level=2) # about 7MB to download

# Subset the SMALLEST state - only 75 municipalities
br_sub <- br[br$NAME_1=="Sergipe" & !is.na(br$NAME_1),]
plot(br_sub)

# Now let's extract the data for each municipality
beginCluster()
e <- extract(b, br_sub, fun=mean, df=T)
endCluster()


Even using the smallest state possible, this example takes about 20 minutes to 
run on a dual-core 2.5GHz Macbook Pro using four threads. As a reference, there 
are states with over 850 municipalities. And remember, in total there are over 
5500 municipalities I need to extract the data for.

I have the feeling that my code is not very optimized.

Has anybody ever dealt with this amount of data in this kind of raster 
operation? Is there any fastest way to achieve my expected result?
Thanks in advance,
-- Thiago V. dos Santos

PhD student
Land and Atmospheric Science

University of Minnesota

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] different models

2016-04-06 Thread Virginia Morera Pujol 
Hi all,

This might be a very dumb question that shows I have very little idea of
what I am talking about, but I'll risk it:

What is the difference between fitting a model using these 3 different
syntaxes?

1/ fit1 <- ppm(ppp, ~covariate),
2/ fit2 <- ppm(ppp, ~x+y+Z, covariates=list(Z=covariate))
3/ fit3 <- ppm(ppp, ~x+y+covariate)

where ppp is my point pattern and "covariate" is a pixel image? I realise
the outputs of 2 and 3 are the same and different to that of 1, so I guess
the question really is

a/ Is there any difference, practical or in the actual computations of the
model, between using 2 and 3?
b/ What is the difference between (2&3) and 1?

Thanks a lot,

Virginia Morera
PhD Student
Department of Animal Biology
University of Barcelona

Aquest correu electrònic i els annexos poden contenir informació
confidencial o protegida legalment i està adreçat exclusivament a la
persona o entitat destinatària. Si no sou  el destinatari final o la
persona encarregada de rebre’l, no esteu autoritzat a llegir-lo,
retenir-lo, modificar-lo, distribuir-lo, copiar-lo ni a revelar-ne el
contingut. Si heu rebut aquest correu electrònic per error, us preguem que
n’informeu al remitent i que elimineu del sistema el missatge i el material
annex que pugui contenir. Gràcies per la vostra col·laboració.

Este correo electrónico y sus anexos pueden contener información
confidencial o legalmente protegida y está exclusivamente dirigido a la
persona o entidad destinataria. Si usted no es el destinatario final o la
persona encargada de recibirlo, no está autorizado a leerlo, retenerlo,
modificarlo, distribuirlo, copiarlo ni a revelar su contenido. Si ha
recibido este mensaje electrónico por error, le rogamos que informe al
remitente y elimine del sistema el mensaje y el material anexo que pueda
contener. Gracias por su colaboración.

This email message and any documents attached to it may contain
confidential or legally protected material and are intended solely for the
use of the individual or organization to whom they are addressed. We remind
you that if you are not the intended recipient of this email message or the
person responsible for processing it, then you are not authorized to read,
save, modify, send, copy or disclose any of its contents. If you have
received this email message by mistake, we kindly ask you to inform the
sender of this and to eliminate both the message and any attachments it
carries from your account.Thank you for your collaboration.

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo