[R-sig-Geo] 8. Convert from spatialpointsdataframe to spatialpoints
Alexandre dos Santos
alexandresantosbr at yahoo.com.br
Fri Nov 5 12:58:56 CET 2010
Dr.Manuel Spimnola,
For create SpatialPoints in shapfile format with
SpatialPointsDataFrame, I use de arguments below,
require(sp)
pontos<-data.frame(p_internos)
head(pontos)
class(pontos)
dim(pontos)
coordinates(pontos) <- c("N_S","L_O")
class(pontos)
dim(pontos)
getSlots("SpatialPointsDataFrame")
slot(pontos, "bbox")
require(maptools)
writePointsShape(pontos,"AMOSTRAGEM", factor2char = TRUE, max_nchar=254)
readShapePoints("AMOSTRAGEM", proj4string = CRS(as.character(NA)),verbose =
FALSE,repair=FALSE)
I will expected help you,
Best regards,
Alexandre dos Santos
Engenheiro Florestal, Msc.
Doutorando em Entomologia
Manejo Integrado de Pragas Florestais
Universidade Federal de Lavras
Departamento de Entomologia
Caixa Postal 3037
37200-000 - Lavras - MG - Brasil
Telefones: +55(35)3829-5122/+55(35)9182-0709
-----Mensagem original-----
De: r-sig-geo-bounces at stat.math.ethz.ch
[mailto:r-sig-geo-bounces at stat.math.ethz.ch] Em nome de
r-sig-geo-request at stat.math.ethz.ch
Enviada em: sexta-feira, 5 de novembro de 2010 09:00
Para: r-sig-geo at stat.math.ethz.ch
Assunto: R-sig-Geo Digest, Vol 87, Issue 5
Send R-sig-Geo mailing list submissions to
r-sig-geo at stat.math.ethz.ch
To subscribe or unsubscribe via the World Wide Web, visit
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
or, via email, send a message with subject or body 'help' to
r-sig-geo-request at stat.math.ethz.ch
You can reach the person managing the list at
r-sig-geo-owner at stat.math.ethz.ch
When replying, please edit your Subject line so it is more specific
than "Re: Contents of R-sig-Geo digest..."
Today's Topics:
1. how to dismiss points outside of my area of interest in a
raster (Manuel Sp?nola)
2. [R-sig-geo]: Plot spatial data with spplot and maps (pelt)
3. Re: [R-sig-geo]: Plot spatial data with spplot and maps
(Karl Ove Hufthammer)
4. raster package and predict for GBM models (Michael ODonnell)
5. Re: how to dismiss points outside of my area of interest in a
raster (Robert J. Hijmans)
6. Re: raster package and predict for GBM models (Robert J. Hijmans)
7. Re: how to dismiss points outside of my area of interest in a
raster (Manuel Sp?nola)
8. Convert from spatialpointsdataframe to spatialpoints
(Manuel Sp?nola)
9. Re: Convert from spatialpointsdataframe to spatialpoints
(Robert J. Hijmans)
10. Re: Convert from spatialpointsdataframe to spatialpoints
(Manuel Sp?nola)
11. Re: Convert from spatialpointsdataframe to spatialpoints
(Robert J. Hijmans)
12. plotting time series for particular months! (govindas at msu.edu)
13. Re: raster package and predict for GBM models (Michael ODonnell)
14. Re: plotting time series for particular months! (steven mosher)
15. Re: Convert from spatialpointsdataframe to spatialpoints
(Manuel Sp?nola)
16. Re: plotting time series for particular months! (govindas at msu.edu)
17. Re: plotting time series for particular months! (steven mosher)
18. rgdal, multiband CCITT compression (Oliver Soong)
19. Re: rgdal, multiband CCITT compression (Michael Sumner)
20. Is there a function for counting polygons within a
multi-polygon SpatialPolygonDataFrame class (Craig Mundy)
21. Re: rgdal, multiband CCITT compression (Oliver Soong)
22. Re: rgdal, multiband CCITT compression (Michael Sumner)
23. Re: rgdal, multiband CCITT compression (Michael Sumner)
24. Re: Is there a function for counting polygons within a
multi-polygon SpatialPolygonDataFrame class (Michael Sumner)
25. Re: rgdal, multiband CCITT compression (Oliver Soong)
26. Defining Holes Inside a Polygon (Rodrigo Aluizio)
27. Re: Defining Holes Inside a Polygon (Roger Bivand)
----------------------------------------------------------------------
Message: 1
Date: Thu, 04 Nov 2010 06:43:58 -0600
From: Manuel Sp?nola <mspinola10 at gmail.com>
To: r-sig-geo at stat.math.ethz.ch
Subject: [R-sig-Geo] how to dismiss points outside of my area of
interest in a raster
Message-ID: <4CD2AA8E.7090804 at gmail.com>
Content-Type: text/plain
Dear list members,
How can eliminate or dismiss points that do not "touch" my raster file?
This mean those points that are outside of the raster area.
When I have my plot of the rester and I draw the points:
points (villages)
some of them are outside of the raster not because an error but just
because the point file is bigger in extent than the raster.
Thank you very much in advance.
Best,
Manuel
--
*Manuel Spmnola, Ph.D.*
Instituto Internacional en Conservacisn y Manejo de Vida Silvestre
Universidad Nacional
Apartado 1350-3000
Heredia
COSTA RICA
mspinola at una.ac.cr
mspinola10 at gmail.com
Telifono: (506) 2277-3598
Fax: (506) 2237-7036
Personal website: Lobito de rmo
<https://sites.google.com/site/lobitoderio/>
Institutional website: ICOMVIS <http://www.icomvis.una.ac.cr/>
[[alternative HTML version deleted]]
------------------------------
Message: 2
Date: Thu, 04 Nov 2010 13:24:46 +0000
From: pelt <pelt at knmi.nl>
To: r-sig-geo at stat.math.ethz.ch
Subject: [R-sig-Geo] [R-sig-geo]: Plot spatial data with spplot and
maps
Message-ID: <4CD2B41E.9050807 at knmi.nl>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Hi all,
I have made a plot with ssplot, using a SpatialPointsDataFrame. The
content is quite simple, as I have 9 grid points with lon/lat
coordinates and 9 values attached to these coordinates. They are in a
square area of 3 by 3 gridboxes.
I would like to lay a map from maps() over these values, but when I try
this, the grids of the maps do not overlap with the grids I have
already created with spplot(). I would like shaded gridboxes, the boxes
should be filled with colour. Can anyone help me with this problem?
[[elided Yahoo spam]]
This is part of my code:
lonlim <-c(4.5,12)
latlim <-c(46.5,52.5)
map(xlim=lonlim, ylim=latlim,col='white')
gt=GridTopology(cellcentre.offset=c(5.75,47.5),cellsize=c(2.5,2),cells.dim=c
(3,3))
grd=SpatialGrid(gt,proj4string = CRS(as.character(NA)))
gridparameters(grd)
gd<-as.data.frame(grd,"SpatialGrid")
h<-season_djf # which are 9 values
h<-as.data.frame(h)
djf.att=SpatialPointsDataFrame(gd,h)
gridded(djf.att)=TRUE
spplot(djf.att, col.regions = colorRampPalette(c("red","orange",
"yellow", "lightblue", "blue",
"purple")), xlab="Longitude(?)",
ylab="Latitude(?)",add=T,scales=list(draw=T),main="Seasonal change",
font.main=4)
map('rivers',add=T)
map('worldHires',add=T)
Regards,
Saskia van Pelt
------------------------------
Message: 3
Date: Thu, 04 Nov 2010 15:32:25 +0100
From: Karl Ove Hufthammer <karl at huftis.org>
To: r-sig-geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] [R-sig-geo]: Plot spatial data with spplot
and maps
Message-ID: <iaug5p$fpe$1 at dough.gmane.org>
Content-Type: text/plain; charset="UTF-8"
pelt wrote:
> I have made a plot with ssplot, using a SpatialPointsDataFrame. The
> content is quite simple, as I have 9 grid points with lon/lat
> coordinates and 9 values attached to these coordinates. They are in a
> square area of 3 by 3 gridboxes.
> I would like to lay a map from maps() over these values, but when I try
> this, the grids of the maps do not overlap with the grids I have
> already created with spplot().
You?re trying to mix ?base? graphics and ?grid? graphics, which doesn?t
work.
?spplot? use ?lattice? graphics (based on ?grid?), so you?ll have to modify
its ?panel? function. Here?s a very simple example.
BTW, I recommend reading the ?lattice? book to understand how to use ?panel?
functions: http://lmdvr.r-forge.r-project.org/figures/figures.html
library(sp)
library(lattice)
library(maps)
library(maptools)
n=100
d=data.frame(x=rnorm(n,0,20),y=rnorm(n,0,20),z=1:n)
xl=range(d$x)
yl=range(d$y)
coordinates(d)=~x+y
spplot(d, panel=function(...) {
panel.fill("lightblue3")
wmap=map("world", fill=TRUE, plot=FALSE,
xlim=xl, ylim=yl, resolution=0)
wmap[c("x", "y")]=mapthin(wmap[c("x", "y")],1)
wmap.sp=map2SpatialPolygons(wmap, wmap$names)
sp.polygons(wmap.sp,fill="khaki", col=NA,lwd=2)
panel.pointsplot(...)
}, xlim=xl, ylim=yl )
The reason for the slow drawing is that each polygon is plotted separately,
and the map polygon also contains areas outside ?xl? and ?yl?. (This is easy
enough to fix by merging and clipping the polygons.)
--
Karl Ove Hufthammer
------------------------------
Message: 4
Date: Thu, 4 Nov 2010 07:34:28 -0700 (PDT)
From: Michael ODonnell <odonnems at yahoo.com>
To: r-sig-geo at stat.math.ethz.ch
Subject: [R-sig-Geo] raster package and predict for GBM models
Message-ID: <519902.86587.qm at web58007.mail.re3.yahoo.com>
Content-Type: text/plain; charset=us-ascii
Hello,
I have little experience using the raster package and I am currently trying
to develop a predicted surface based on a gradient boosted model (gbm
package).
The basic steps I am using include the following:
# Creating a raster stack of independent variables
USHGT <- raster("ushgt.asc", package="raster")
... (other variables omitted for post)
ras_stk <- stack(USHGT, HYELB, T1TO5LE, URBD_D, SG_UNIT, MINTEMP, ANNPCP,
TC_B04, TC_B07)
# Create the prediction model
... I am using a separate script that creates the gbm object and due to the
length I will not include. The outcome is a gbm object.
.. Now I use the following, which fails:
r <- predict(ras_stk, myData, filename=OutFile, na.rm=TRUE, overwrite=TRUE,
n.trees=myData$n.trees, type="response", progress="window")
Error:
Error in model.frame.default(terms(reformulate(object$var.names)), newdata,
:object is not a matrix
The raster package manual states that the predict function will work on any
model that a predict method has been implemented.
Does anyone have experience with this or possibly an idea on how I might
solve the issue.
Thank you,
mike
------------------------------
Message: 5
Date: Thu, 4 Nov 2010 09:46:15 -0700
From: "Robert J. Hijmans" <r.hijmans at gmail.com>
To: Manuel Sp?nola <mspinola10 at gmail.com>
Cc: r-sig-geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] how to dismiss points outside of my area of
interest in a raster
Message-ID:
<AANLkTi=7LEC_wLjJSdrQLpkhaC7yRCnEO+sCr2bO9nxF at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
Manuel,
library(raster)
# Assuming you have matrix of points (or SpatialPoints*)
pt <- (cbind(runif(100), runif(100)) -0.5) * 500
spt <- SpatialPoints(pt)
# and a RasterLayer
r <- raster()
# you can do
e <- polygonFromExtent(r)
o <- overlay(e, spt)
o <- which(!is.na(o))
pts <- pt[o,]
# show that it worked
plot(pt)
r[] = 1
plot(r, add=TRUE, legend=FALSE)
points(pt, col='blue')
points(pts, col='red', pch=20)
Robert
On Thu, Nov 4, 2010 at 5:43 AM, Manuel Sp?nola <mspinola10 at gmail.com> wrote:
> Dear list members,
>
> How can eliminate or dismiss points that do not "touch" my raster file?
> This mean those points that are outside of the raster area.
> When I have my plot of the rester and I draw the points:
>
> points (villages)
>
> some of them are outside of the raster not because an error but just
> because the point file is bigger in extent than the raster.
>
> Thank you very much in advance.
> Best,
>
> Manuel
>
>
> --
> *Manuel Sp?nola, Ph.D.*
> Instituto Internacional en Conservaci?n y Manejo de Vida Silvestre
> Universidad Nacional
> Apartado 1350-3000
> Heredia
> COSTA RICA
> mspinola at una.ac.cr
> mspinola10 at gmail.com
> Tel?fono: (506) 2277-3598
> Fax: (506) 2237-7036
> 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
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
------------------------------
Message: 6
Date: Thu, 4 Nov 2010 09:58:37 -0700
From: "Robert J. Hijmans" <r.hijmans at gmail.com>
To: Michael ODonnell <odonnems at yahoo.com>
Cc: r-sig-geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] raster package and predict for GBM models
Message-ID:
<AANLkTi=7TBsyUVWcOLiGqEi_tTRwDNSDcjaRVXODKV-t at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
Hi Mike,
I cannot know from what you send, but the second argument should be
your gbm object, and from your name "myData" I am guessing it is not
(unless it really is referring to youModel).
One of the two vignettes of the 'dismo' package is about 'gbm'
(boosted regression trees) and 'raster' . So in principle it should
work.
Robert
On Thu, Nov 4, 2010 at 7:34 AM, Michael ODonnell <odonnems at yahoo.com> wrote:
> Hello,
>
> I have little experience using the raster package and I am currently
trying to develop a predicted surface based on a gradient boosted model (gbm
package).
>
> The basic steps I am using include the following:
>
> # Creating a raster stack of independent variables
> USHGT <- raster("ushgt.asc", package="raster")
> ... (other variables omitted for post)
> ras_stk <- stack(USHGT, HYELB, T1TO5LE, URBD_D, SG_UNIT, MINTEMP, ANNPCP,
TC_B04, TC_B07)
>
> # Create the prediction model
> ... I am using a separate script that creates the gbm object and due to
the length I will not include. The outcome is a gbm object.
>
> .. Now I use the following, which fails:
> r <- predict(ras_stk, myData, filename=OutFile, na.rm=TRUE,
overwrite=TRUE, n.trees=myData$n.trees, type="response", progress="window")
>
> Error:
> Error in model.frame.default(terms(reformulate(object$var.names)),
newdata, ?:object is not a matrix
>
> The raster package manual states that the predict function will work on
any model that a predict method has been implemented.
>
> Does anyone have experience with this or possibly an idea on how I might
solve the issue.
>
> Thank you,
> mike
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
------------------------------
Message: 7
Date: Thu, 04 Nov 2010 11:21:30 -0600
From: Manuel Sp?nola <mspinola10 at gmail.com>
To: "Robert J. Hijmans" <r.hijmans at gmail.com>
Cc: r-sig-geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] how to dismiss points outside of my area of
interest in a raster
Message-ID: <4CD2EB9A.5040602 at gmail.com>
Content-Type: text/plain
Thank you very much Robert.
Best,
Manuel
On 04/11/2010 10:46 a.m., Robert J. Hijmans wrote:
> Manuel,
>
> library(raster)
> # Assuming you have matrix of points (or SpatialPoints*)
> pt<- (cbind(runif(100), runif(100)) -0.5) * 500
> spt<- SpatialPoints(pt)
> # and a RasterLayer
> r<- raster()
>
> # you can do
>
> e<- polygonFromExtent(r)
> o<- overlay(e, spt)
> o<- which(!is.na(o))
> pts<- pt[o,]
>
>
> # show that it worked
> plot(pt)
> r[] = 1
> plot(r, add=TRUE, legend=FALSE)
> points(pt, col='blue')
> points(pts, col='red', pch=20)
>
>
> Robert
>
> On Thu, Nov 4, 2010 at 5:43 AM, Manuel Spmnola<mspinola10 at gmail.com>
wrote:
>> Dear list members,
>>
>> How can eliminate or dismiss points that do not "touch" my raster file?
>> This mean those points that are outside of the raster area.
>> When I have my plot of the rester and I draw the points:
>>
>> points (villages)
>>
>> some of them are outside of the raster not because an error but just
>> because the point file is bigger in extent than the raster.
>>
>> Thank you very much in advance.
>> Best,
>>
>> Manuel
>>
>>
>> --
>> *Manuel Spmnola, Ph.D.*
>> Instituto Internacional en Conservacisn y Manejo de Vida Silvestre
>> Universidad Nacional
>> Apartado 1350-3000
>> Heredia
>> COSTA RICA
>> mspinola at una.ac.cr
>> mspinola10 at gmail.com
>> Telifono: (506) 2277-3598
>> Fax: (506) 2237-7036
>> Personal website: Lobito de rmo
>> <https://sites.google.com/site/lobitoderio/>
>> Institutional website: ICOMVIS<http://www.icomvis.una.ac.cr/>
>>
>> [[alternative HTML version deleted]]
>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
--
*Manuel Spmnola, Ph.D.*
Instituto Internacional en Conservacisn y Manejo de Vida Silvestre
Universidad Nacional
Apartado 1350-3000
Heredia
COSTA RICA
mspinola at una.ac.cr
mspinola10 at gmail.com
Telifono: (506) 2277-3598
Fax: (506) 2237-7036
Personal website: Lobito de rmo
<https://sites.google.com/site/lobitoderio/>
Institutional website: ICOMVIS <http://www.icomvis.una.ac.cr/>
[[alternative HTML version deleted]]
------------------------------
Message: 8
Date: Thu, 04 Nov 2010 11:25:51 -0600
From: Manuel Sp?nola <mspinola10 at gmail.com>
To: r-sig-geo at stat.math.ethz.ch
Subject: [R-sig-Geo] Convert from spatialpointsdataframe to
spatialpoints
Message-ID: <4CD2EC9F.2050000 at gmail.com>
Content-Type: text/plain
Dear list members,
I am trying to create a raster from points, but I read the points from a
shapfile, so I have the points object as a SpatialPointsDataFrame.
To use the function pointsToRaster I need the points as SpatialPoints.
How to convert a SpatialPointsDataFrame object to a SpatialPoints?
Best,
Manuel
--
*Manuel Spmnola, Ph.D.*
Instituto Internacional en Conservacisn y Manejo de Vida Silvestre
Universidad Nacional
Apartado 1350-3000
Heredia
COSTA RICA
mspinola at una.ac.cr
mspinola10 at gmail.com
Telifono: (506) 2277-3598
Fax: (506) 2237-7036
Personal website: Lobito de rmo
<https://sites.google.com/site/lobitoderio/>
Institutional website: ICOMVIS <http://www.icomvis.una.ac.cr/>
[[alternative HTML version deleted]]
------------------------------
Message: 9
Date: Thu, 4 Nov 2010 10:41:15 -0700
From: "Robert J. Hijmans" <r.hijmans at gmail.com>
To: Manuel Sp?nola <mspinola10 at gmail.com>
Cc: r-sig-geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] Convert from spatialpointsdataframe to
spatialpoints
Message-ID:
<AANLkTikub6cdL3R5WGssPikcOUah0+MedPwqERXPDvme at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
Manuel.
Did you try it without converting? The help file in raster says it
needs a SpatialPoints* object. The '*' is a wild card that in this
case stands for 'DataFrame'.
The answer to your question:
as(x, 'SpatialPoints')
You could find the answer to many of this type of questions by reading
the vignettes of sp and raster, and by searching the R-sig-geo help
archives, and looking at the R spatial wiki.
Robert
On Thu, Nov 4, 2010 at 10:25 AM, Manuel Sp?nola <mspinola10 at gmail.com>
wrote:
> Dear list members,
>
>
> I am trying to create a raster from points, but I read the points from a
> shapfile, so I have the points object as a SpatialPointsDataFrame.
> To use the function pointsToRaster I need the points as SpatialPoints.
>
> How to convert a SpatialPointsDataFrame object to a SpatialPoints?
> Best,
>
> Manuel
>
>
> --
> *Manuel Sp?nola, Ph.D.*
> Instituto Internacional en Conservaci?n y Manejo de Vida Silvestre
> Universidad Nacional
> Apartado 1350-3000
> Heredia
> COSTA RICA
> mspinola at una.ac.cr
> mspinola10 at gmail.com
> Tel?fono: (506) 2277-3598
> Fax: (506) 2237-7036
> 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
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
------------------------------
Message: 10
Date: Thu, 04 Nov 2010 12:02:30 -0600
From: Manuel Sp?nola <mspinola10 at gmail.com>
To: "Robert J. Hijmans" <r.hijmans at gmail.com>
Cc: r-sig-geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] Convert from spatialpointsdataframe to
spatialpoints
Message-ID: <4CD2F536.9050102 at gmail.com>
Content-Type: text/plain
Thank you very much Robert.
I am reading the raster vignette right now.
I got this error:
> r = raster(corredor)
> res(r) = 50
> r = pointsToRaster(poblados, r)
Error en .pointsToMatrix(xy) :
points should be vectors of length 2, matrices with 2 columns, or a
SpatialPoints* object
> poblados at data
x y
1 500577.8 1155145
2 495470.2 1164700
3 489082.0 1152228
4 494559.3 1171489
5 492094.0 1157327
6 490937.9 1155423
7 488748.0 1153979
8 491910.8 1154378
9 485430.7 1148051
10 500060.9 1148109
11 491212.5 1159109
12 487253.9 1151491
13 480244.6 1181636
14 480651.1 1184371
15 477555.2 1181189
16 486482.8 1167946
17 483764.4 1164624
18 477301.6 1192979
19 483737.4 1192048
20 480760.2 1192652
21 476126.6 1165617
22 481182.4 1185941
Best,
Manuel
On 04/11/2010 11:41 a.m., Robert J. Hijmans wrote:
> Manuel.
>
> Did you try it without converting? The help file in raster says it
> needs a SpatialPoints* object. The '*' is a wild card that in this
> case stands for 'DataFrame'.
>
> The answer to your question:
>
> as(x, 'SpatialPoints')
>
> You could find the answer to many of this type of questions by reading
> the vignettes of sp and raster, and by searching the R-sig-geo help
> archives, and looking at the R spatial wiki.
>
> Robert
>
> On Thu, Nov 4, 2010 at 10:25 AM, Manuel Spmnola<mspinola10 at gmail.com>
wrote:
>> Dear list members,
>>
>>
>> I am trying to create a raster from points, but I read the points from a
>> shapfile, so I have the points object as a SpatialPointsDataFrame.
>> To use the function pointsToRaster I need the points as SpatialPoints.
>>
>> How to convert a SpatialPointsDataFrame object to a SpatialPoints?
>> Best,
>>
>> Manuel
>>
>>
>> --
>> *Manuel Spmnola, Ph.D.*
>> Instituto Internacional en Conservacisn y Manejo de Vida Silvestre
>> Universidad Nacional
>> Apartado 1350-3000
>> Heredia
>> COSTA RICA
>> mspinola at una.ac.cr
>> mspinola10 at gmail.com
>> Telifono: (506) 2277-3598
>> Fax: (506) 2237-7036
>> Personal website: Lobito de rmo
>> <https://sites.google.com/site/lobitoderio/>
>> Institutional website: ICOMVIS<http://www.icomvis.una.ac.cr/>
>>
>> [[alternative HTML version deleted]]
>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
--
*Manuel Spmnola, Ph.D.*
Instituto Internacional en Conservacisn y Manejo de Vida Silvestre
Universidad Nacional
Apartado 1350-3000
Heredia
COSTA RICA
mspinola at una.ac.cr
mspinola10 at gmail.com
Telifono: (506) 2277-3598
Fax: (506) 2237-7036
Personal website: Lobito de rmo
<https://sites.google.com/site/lobitoderio/>
Institutional website: ICOMVIS <http://www.icomvis.una.ac.cr/>
[[alternative HTML version deleted]]
------------------------------
Message: 11
Date: Thu, 4 Nov 2010 12:16:47 -0700
From: "Robert J. Hijmans" <r.hijmans at gmail.com>
To: Manuel Sp?nola <mspinola10 at gmail.com>
Cc: r-sig-geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] Convert from spatialpointsdataframe to
spatialpoints
Message-ID:
<AANLkTimW+bPp3v05DfCK1+WbbmbR64c8517YjJXxJizm at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
Manuel,
I think this is wrong:
>> r = pointsToRaster(poblados, r)
Try:
r <- pointsToRaster(r, poblados)
Robert
On Thu, Nov 4, 2010 at 11:02 AM, Manuel Sp?nola <mspinola10 at gmail.com>
wrote:
> Thank you very much Robert.
>
> I am reading the raster vignette right now.
>
> I got this error:
>
>> r = raster(corredor)
>> res(r) = 50
>> r = pointsToRaster(poblados, r)
> Error en .pointsToMatrix(xy) :
> ? points should be vectors of length 2, matrices with 2 columns, or a
> SpatialPoints* object
>> poblados at data
> ????????? x?????? y
> 1? 500577.8 1155145
> 2? 495470.2 1164700
> 3? 489082.0 1152228
> 4? 494559.3 1171489
> 5? 492094.0 1157327
> 6? 490937.9 1155423
> 7? 488748.0 1153979
> 8? 491910.8 1154378
> 9? 485430.7 1148051
> 10 500060.9 1148109
> 11 491212.5 1159109
> 12 487253.9 1151491
> 13 480244.6 1181636
> 14 480651.1 1184371
> 15 477555.2 1181189
> 16 486482.8 1167946
> 17 483764.4 1164624
> 18 477301.6 1192979
> 19 483737.4 1192048
> 20 480760.2 1192652
> 21 476126.6 1165617
> 22 481182.4 1185941
>
>
> Best,
>
> Manuel
>
> On 04/11/2010 11:41 a.m., Robert J. Hijmans wrote:
>
> Manuel.
>
> Did you try it without converting? The help file in raster says it
> needs a SpatialPoints* object. The '*' is a wild card that in this
> case stands for 'DataFrame'.
>
> The answer to your question:
>
> as(x, 'SpatialPoints')
>
> You could find the answer to many of this type of questions by reading
> the vignettes of sp and raster, and by searching the R-sig-geo help
> archives, and looking at the R spatial wiki.
>
> Robert
>
> On Thu, Nov 4, 2010 at 10:25 AM, Manuel Sp?nola <mspinola10 at gmail.com>
> wrote:
>
> Dear list members,
>
>
> I am trying to create a raster from points, but I read the points from a
> shapfile, so I have the points object as a SpatialPointsDataFrame.
> To use the function pointsToRaster I need the points as SpatialPoints.
>
> How to convert a SpatialPointsDataFrame object to a SpatialPoints?
> Best,
>
> Manuel
>
>
> --
> *Manuel Sp?nola, Ph.D.*
> Instituto Internacional en Conservaci?n y Manejo de Vida Silvestre
> Universidad Nacional
> Apartado 1350-3000
> Heredia
> COSTA RICA
> mspinola at una.ac.cr
> mspinola10 at gmail.com
> Tel?fono: (506) 2277-3598
> Fax: (506) 2237-7036
> 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
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
>
>
> --
> Manuel Sp?nola, Ph.D.
> Instituto Internacional en Conservaci?n y Manejo de Vida Silvestre
> Universidad Nacional
> Apartado 1350-3000
> Heredia
> COSTA RICA
> mspinola at una.ac.cr
> mspinola10 at gmail.com
> Tel?fono: (506) 2277-3598
> Fax: (506) 2237-7036
> Personal website: Lobito de r?o
> Institutional website: ICOMVIS
------------------------------
Message: 12
Date: Thu, 04 Nov 2010 15:18:53 -0400
From: govindas at msu.edu
To: r-sig-geo at stat.math.ethz.ch
[[elided Yahoo spam]]
Message-ID: <20101104151853.93076j1dwbv9pwrh at mail.msu.edu>
Content-Type: text/plain
Hi all,
I have a matrix as given below...
year month day prec
1 1980 10 1 13.4
2 1980 10 2 0.0
3 1980 10 3 0.0
4 1980 10 4 0.0
5 1980 10 5 0.0
6 1980 10 6 0.0
7 1980 10 7 0.0
8 1980 10 8 0.0
9 1980 10 9 0.0
10 1980 10 10 0.0
11 1980 10 11 7.4
12 1980 10 12 5.4
13 1980 10 13 7.2
14 1980 10 14 0.0
15 1980 10 15 0.0
16 1980 10 16 0.0
17 1980 10 17 41.2
18 1980 10 18 0.0
19 1980 10 19 0.0
20 1980 10 20 0.0
21 1980 10 21 0.0
22 1980 10 22 0.0
23 1980 10 23 0.0
24 1980 10 24 0.0
25 1980 10 25 0.0
26 1980 10 26 0.0
27 1980 10 27 2.0
28 1980 10 28 0.0
29 1980 10 29 0.0
30 1980 10 30 0.0
31 1980 10 31 0.0
32 1980 11 1 0.0
33 1980 11 2 0.0
34 1980 11 3 0.0
35 1980 11 4 0.0
36 1980 11 5 12.4
the precipitation values extend from 1980 to 2005, but only for october,
november and december. I would like to plot just these 3 months for the
given time period (1980 - 2005). Is there a way to get these values in the
x-axis. i.e. i need a plot with its axis reading "1980, 1981, 1982 .. 2005"
(even if the months are specified then it should be of more use). for now,
i get an axis like .. "0, 500, ... 2000" (the plot is giving the index
values).
[[elided Yahoo spam]]
ts.chn <- ts(chn.arr[1:2386,4], start=c(1980, 10), end=c(2005, 12),
freq=365)
plot(ts.chn)
--
Regards,
Mahalakshmi
Graduate Student
#20, Department of Geography
Michigan State University
East Lansing, MI 48824
[[alternative HTML version deleted]]
------------------------------
Message: 13
Date: Thu, 4 Nov 2010 12:26:45 -0700 (PDT)
From: Michael ODonnell <odonnems at yahoo.com>
To: "Robert J. Hijmans" <r.hijmans at gmail.com>
Cc: r-sig-geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] raster package and predict for GBM models
Message-ID: <468619.42773.qm at web58003.mail.re3.yahoo.com>
Content-Type: text/plain; charset=iso-8859-1
Robert, sorry. Yes, myData is my gbm object.
Also, I figured out what I had done wrong. The raster stack I created
specified a data set name that was lowercase and the the gbm model was
specifying upper case for the same raster layer. Once I use the correct
name, it worked.
Thank you for the response.
Mike
--- On Thu, 11/4/10, Robert J. Hijmans <r.hijmans at gmail.com> wrote:
> From: Robert J. Hijmans <r.hijmans at gmail.com>
> Subject: Re: [R-sig-Geo] raster package and predict for GBM models
> To: "Michael ODonnell" <odonnems at yahoo.com>
> Cc: r-sig-geo at stat.math.ethz.ch
> Date: Thursday, November 4, 2010, 10:58 AM
> Hi Mike,
>
> I cannot know from what you send, but the second argument
> should be
> your gbm object, and from your name "myData" I am guessing
> it is not
> (unless it really is referring to youModel).
> One of the two vignettes of the 'dismo' package is about
> 'gbm'
> (boosted regression trees) and 'raster' . So in principle
> it should
> work.
>
> Robert
>
>
> On Thu, Nov 4, 2010 at 7:34 AM, Michael ODonnell <odonnems at yahoo.com>
> wrote:
> > Hello,
> >
> > I have little experience using the raster package and
> I am currently trying to develop a predicted surface based
> on a gradient boosted model (gbm package).
> >
> > The basic steps I am using include the following:
> >
> > # Creating a raster stack of independent variables
> > USHGT <- raster("ushgt.asc", package="raster")
> > ... (other variables omitted for post)
> > ras_stk <- stack(USHGT, HYELB, T1TO5LE, URBD_D,
> SG_UNIT, MINTEMP, ANNPCP, TC_B04, TC_B07)
> >
> > # Create the prediction model
> > ... I am using a separate script that creates the gbm
> object and due to the length I will not include. The outcome
> is a gbm object.
> >
> > .. Now I use the following, which fails:
> > r <- predict(ras_stk, myData, filename=OutFile,
> na.rm=TRUE, overwrite=TRUE, n.trees=myData$n.trees,
> type="response", progress="window")
> >
> > Error:
> > Error in
> model.frame.default(terms(reformulate(object$var.names)),
> newdata, ?:object is not a matrix
> >
> > The raster package manual states that the predict
> function will work on any model that a predict method has
> been implemented.
> >
> > Does anyone have experience with this or possibly an
> idea on how I might solve the issue.
> >
> > Thank you,
> > mike
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo at stat.math.ethz.ch
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
------------------------------
Message: 14
Date: Thu, 4 Nov 2010 12:33:04 -0700
From: steven mosher <moshersteven at gmail.com>
To: govindas at msu.edu
Cc: r-sig-geo at stat.math.ethz.ch
[[elided Yahoo spam]]
Message-ID:
<AANLkTinWf+hi+y_8bpZ8znk5-WNWyTY6Y3JSAQOvD01u at mail.gmail.com>
Content-Type: text/plain
I would look at using the package 'zoo' your time series data would be a
zoo object
and then your date data would be coerced as well,
zoo, is very slick.
On Thu, Nov 4, 2010 at 12:18 PM, <govindas at msu.edu> wrote:
>
>
> Hi all,
>
> I have a matrix as given below...
>
> year month day prec
> 1 1980 10 1 13.4
> 2 1980 10 2 0.0
> 3 1980 10 3 0.0
> 4 1980 10 4 0.0
> 5 1980 10 5 0.0
> 6 1980 10 6 0.0
> 7 1980 10 7 0.0
> 8 1980 10 8 0.0
> 9 1980 10 9 0.0
> 10 1980 10 10 0.0
> 11 1980 10 11 7.4
> 12 1980 10 12 5.4
> 13 1980 10 13 7.2
> 14 1980 10 14 0.0
> 15 1980 10 15 0.0
> 16 1980 10 16 0.0
> 17 1980 10 17 41.2
> 18 1980 10 18 0.0
> 19 1980 10 19 0.0
> 20 1980 10 20 0.0
> 21 1980 10 21 0.0
> 22 1980 10 22 0.0
> 23 1980 10 23 0.0
> 24 1980 10 24 0.0
> 25 1980 10 25 0.0
> 26 1980 10 26 0.0
> 27 1980 10 27 2.0
> 28 1980 10 28 0.0
> 29 1980 10 29 0.0
> 30 1980 10 30 0.0
> 31 1980 10 31 0.0
> 32 1980 11 1 0.0
> 33 1980 11 2 0.0
> 34 1980 11 3 0.0
> 35 1980 11 4 0.0
> 36 1980 11 5 12.4
>
> the precipitation values extend from 1980 to 2005, but only for october,
> november and december. I would like to plot just these 3 months for the
> given time period (1980 - 2005). Is there a way to get these values in
the
> x-axis. i.e. i need a plot with its axis reading "1980, 1981, 1982 ..
2005"
> (even if the months are specified then it should be of more use). for
now,
> i get an axis like .. "0, 500, ... 2000" (the plot is giving the index
> values).
>
[[elided Yahoo spam]]
> ts.chn <- ts(chn.arr[1:2386,4], start=c(1980, 10), end=c(2005, 12),
> freq=365)
> plot(ts.chn)
>
> --
> Regards,
> Mahalakshmi
> Graduate Student
> #20, Department of Geography
> Michigan State University
> East Lansing, MI 48824
> [[alternative HTML version deleted]]
>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
[[alternative HTML version deleted]]
------------------------------
Message: 15
Date: Thu, 04 Nov 2010 14:07:27 -0600
From: Manuel Sp?nola <mspinola10 at gmail.com>
To: "Robert J. Hijmans" <r.hijmans at gmail.com>
Cc: r-sig-geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] Convert from spatialpointsdataframe to
spatialpoints
Message-ID: <4CD3127F.9040309 at gmail.com>
Content-Type: text/plain
Thank you very much Robert,
It works but it looks empty.
I don't have values, is only xy values (the locations for villages).
I did:
r = pointsToRaster(r, poblados, values=1)
but still looks empty.
Best,
Manuel
On 04/11/2010 01:16 p.m., Robert J. Hijmans wrote:
> Manuel,
>
> I think this is wrong:
>>> r = pointsToRaster(poblados, r)
> Try:
> r<- pointsToRaster(r, poblados)
>
>
> Robert
>
> On Thu, Nov 4, 2010 at 11:02 AM, Manuel Spmnola<mspinola10 at gmail.com>
wrote:
>> Thank you very much Robert.
>>
>> I am reading the raster vignette right now.
>>
>> I got this error:
>>
>>> r = raster(corredor)
>>> res(r) = 50
>>> r = pointsToRaster(poblados, r)
>> Error en .pointsToMatrix(xy) :
>> points should be vectors of length 2, matrices with 2 columns, or a
>> SpatialPoints* object
>>> poblados at data
>> x y
>> 1 500577.8 1155145
>> 2 495470.2 1164700
>> 3 489082.0 1152228
>> 4 494559.3 1171489
>> 5 492094.0 1157327
>> 6 490937.9 1155423
>> 7 488748.0 1153979
>> 8 491910.8 1154378
>> 9 485430.7 1148051
>> 10 500060.9 1148109
>> 11 491212.5 1159109
>> 12 487253.9 1151491
>> 13 480244.6 1181636
>> 14 480651.1 1184371
>> 15 477555.2 1181189
>> 16 486482.8 1167946
>> 17 483764.4 1164624
>> 18 477301.6 1192979
>> 19 483737.4 1192048
>> 20 480760.2 1192652
>> 21 476126.6 1165617
>> 22 481182.4 1185941
>>
>>
>> Best,
>>
>> Manuel
>>
>> On 04/11/2010 11:41 a.m., Robert J. Hijmans wrote:
>>
>> Manuel.
>>
>> Did you try it without converting? The help file in raster says it
>> needs a SpatialPoints* object. The '*' is a wild card that in this
>> case stands for 'DataFrame'.
>>
>> The answer to your question:
>>
>> as(x, 'SpatialPoints')
>>
>> You could find the answer to many of this type of questions by reading
>> the vignettes of sp and raster, and by searching the R-sig-geo help
>> archives, and looking at the R spatial wiki.
>>
>> Robert
>>
>> On Thu, Nov 4, 2010 at 10:25 AM, Manuel Spmnola<mspinola10 at gmail.com>
>> wrote:
>>
>> Dear list members,
>>
>>
>> I am trying to create a raster from points, but I read the points from a
>> shapfile, so I have the points object as a SpatialPointsDataFrame.
>> To use the function pointsToRaster I need the points as SpatialPoints.
>>
>> How to convert a SpatialPointsDataFrame object to a SpatialPoints?
>> Best,
>>
>> Manuel
>>
>>
>> --
>> *Manuel Spmnola, Ph.D.*
>> Instituto Internacional en Conservacisn y Manejo de Vida Silvestre
>> Universidad Nacional
>> Apartado 1350-3000
>> Heredia
>> COSTA RICA
>> mspinola at una.ac.cr
>> mspinola10 at gmail.com
>> Telifono: (506) 2277-3598
>> Fax: (506) 2237-7036
>> Personal website: Lobito de rmo
>> <https://sites.google.com/site/lobitoderio/>
>> Institutional website: ICOMVIS<http://www.icomvis.una.ac.cr/>
>>
>> [[alternative HTML version deleted]]
>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
>>
>>
>> --
>> Manuel Spmnola, Ph.D.
>> Instituto Internacional en Conservacisn y Manejo de Vida Silvestre
>> Universidad Nacional
>> Apartado 1350-3000
>> Heredia
>> COSTA RICA
>> mspinola at una.ac.cr
>> mspinola10 at gmail.com
>> Telifono: (506) 2277-3598
>> Fax: (506) 2237-7036
>> Personal website: Lobito de rmo
>> Institutional website: ICOMVIS
--
*Manuel Spmnola, Ph.D.*
Instituto Internacional en Conservacisn y Manejo de Vida Silvestre
Universidad Nacional
Apartado 1350-3000
Heredia
COSTA RICA
mspinola at una.ac.cr
mspinola10 at gmail.com
Telifono: (506) 2277-3598
Fax: (506) 2237-7036
Personal website: Lobito de rmo
<https://sites.google.com/site/lobitoderio/>
Institutional website: ICOMVIS <http://www.icomvis.una.ac.cr/>
[[alternative HTML version deleted]]
------------------------------
Message: 16
Date: Thu, 04 Nov 2010 17:36:22 -0400
From: govindas at msu.edu
To: steven mosher <moshersteven at gmail.com>
Cc: r-sig-geo at stat.math.ethz.ch
[[elided Yahoo spam]]
Message-ID: <20101104173622.327465q203o894k6 at mail.msu.edu>
Content-Type: text/plain; charset="iso-8859-1"
thanks .. and it did work! but i see some jagged lines in between the years,
typically like the one? joining dec of previous year and oct of the
successive year. Is this how it would look like? can it be made better?
[[elided Yahoo spam]]
heres the code..!
x.date <- as.Date(paste(chn.arr[,1], chn.arr[,2], chn.arr[,3], sep="-"))
ts.chn <- ts(chn.arr[1:2386,4])
z.chn <- zoo(ts.chn, x.date)
plot(z.chn, main = "Rainfall variations", xlab = "Year", ylab="Daily
rainfall (mm)")
--
Regards,
Mahalakshmi
Graduate Student
#20, Department of Geography
Michigan State University
East Lansing, MI 48824 Quoting steven mosher <moshersteven at gmail.com>:
> I would look at using? the package 'zoo' your time series data would be a
> zoo object
> and then your date data? would be coerced as well,
>
> zoo, is very slick.
>
> On Thu, Nov 4, 2010 at 12:18 PM, <govindas at msu.edu> wrote:
>
>>
>>
>> Hi all,
>>
>> I have a matrix as given below...
>>
>> year month day prec
>> 1? ?1980? ? 10? ?1 13.4
>> 2? ?1980? ? 10? ?2? 0.0
>> 3? ?1980? ? 10? ?3? 0.0
>> 4? ?1980? ? 10? ?4? 0.0
>> 5? ?1980? ? 10? ?5? 0.0
>> 6? ?1980? ? 10? ?6? 0.0
>> 7? ?1980? ? 10? ?7? 0.0
>> 8? ?1980? ? 10? ?8? 0.0
>> 9? ?1980? ? 10? ?9? 0.0
>> 10? 1980? ? 10? 10? 0.0
>> 11? 1980? ? 10? 11? 7.4
>> 12? 1980? ? 10? 12? 5.4
>> 13? 1980? ? 10? 13? 7.2
>> 14? 1980? ? 10? 14? 0.0
>> 15? 1980? ? 10? 15? 0.0
>> 16? 1980? ? 10? 16? 0.0
>> 17? 1980? ? 10? 17 41.2
>> 18? 1980? ? 10? 18? 0.0
>> 19? 1980? ? 10? 19? 0.0
>> 20? 1980? ? 10? 20? 0.0
>> 21? 1980? ? 10? 21? 0.0
>> 22? 1980? ? 10? 22? 0.0
>> 23? 1980? ? 10? 23? 0.0
>> 24? 1980? ? 10? 24? 0.0
>> 25? 1980? ? 10? 25? 0.0
>> 26? 1980? ? 10? 26? 0.0
>> 27? 1980? ? 10? 27? 2.0
>> 28? 1980? ? 10? 28? 0.0
>> 29? 1980? ? 10? 29? 0.0
>> 30? 1980? ? 10? 30? 0.0
>> 31? 1980? ? 10? 31? 0.0
>> 32? 1980? ? 11? ?1? 0.0
>> 33? 1980? ? 11? ?2? 0.0
>> 34? 1980? ? 11? ?3? 0.0
>> 35? 1980? ? 11? ?4? 0.0
>> 36? 1980? ? 11? ?5 12.4
>>
>> the precipitation values extend from 1980 to 2005, but only for october,
>> november and december. I would? like to plot just these 3 months for the
>> given time period (1980 -? 2005). Is there a way to get these values in
the
>> x-axis. i.e. i need a? plot with its axis reading "1980, 1981, 1982 ..
2005"
>> (even if the? months are specified then it should be of more use). for
now,
>> i get an? axis like .. "0, 500, ... 2000" (the plot is giving the index
>> values).
>>
[[elided Yahoo spam]]
>> ts.chn <- ts(chn.arr[1:2386,4], start=c(1980, 10), end=c(2005, 12),
>> freq=365)
>> plot(ts.chn)
>>
>> --
>> Regards,
>> Mahalakshmi
>> Graduate Student
>> #20, Department of Geography
>> Michigan State University
>> East Lansing, MI 48824
>>? ? ? ? [[alternative HTML version deleted]]
>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL:
<https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20101104/f6396585/atta
chment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: RF_zoo.jpeg
Type: image/jpeg
Size: 79123 bytes
Desc: not available
URL:
<https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20101104/f6396585/atta
chment-0001.jpeg>
------------------------------
Message: 17
Date: Thu, 4 Nov 2010 15:25:05 -0700
From: steven mosher <moshersteven at gmail.com>
To: govindas at msu.edu
Cc: r-sig-geo at stat.math.ethz.ch
[[elided Yahoo spam]]
Message-ID:
<AANLkTimgiRAHchqE361VRRk8u6SeZYeddhWobM+Lgabw at mail.gmail.com>
Content-Type: text/plain
Hmm,
Suggests you might have some missing data as opposed to having zeros for
those days.
something about your number of days for oct/nov/dec from 1980-2005 doesnt
seem right either
so double check that.
On Thu, Nov 4, 2010 at 2:36 PM, <govindas at msu.edu> wrote:
> thanks .. and it did work! but i see some jagged lines in between the
> years, typically like the one joining dec of previous year and oct of the
> successive year. Is this how it would look like? can it be made better?
>
[[elided Yahoo spam]]
>
> heres the code..!
>
> x.date <- as.Date(paste(chn.arr[,1], chn.arr[,2], chn.arr[,3], sep="-"))
> ts.chn <- ts(chn.arr[1:2386,4])
> z.chn <- zoo(ts.chn, x.date)
> plot(z.chn, main = "Rainfall variations", xlab = "Year", ylab="Daily
> rainfall (mm)")
>
> --
> Regards,
> Mahalakshmi
> Graduate Student
> #20, Department of Geography
> Michigan State University
> East Lansing, MI 48824
> Quoting steven mosher <moshersteven at gmail.com>:
>
> > I would look at using the package 'zoo' your time series data would be
a
> > zoo object
> > and then your date data would be coerced as well,
> >
> > zoo, is very slick.
> >
> > On Thu, Nov 4, 2010 at 12:18 PM, <govindas at msu.edu> wrote:
> >
> >>
> >>
> >> Hi all,
> >>
> >> I have a matrix as given below...
> >>
> >> year month day prec
> >> 1 1980 10 1 13.4
> >> 2 1980 10 2 0.0
> >> 3 1980 10 3 0.0
> >> 4 1980 10 4 0.0
> >> 5 1980 10 5 0.0
> >> 6 1980 10 6 0.0
> >> 7 1980 10 7 0.0
> >> 8 1980 10 8 0.0
> >> 9 1980 10 9 0.0
> >> 10 1980 10 10 0.0
> >> 11 1980 10 11 7.4
> >> 12 1980 10 12 5.4
> >> 13 1980 10 13 7.2
> >> 14 1980 10 14 0.0
> >> 15 1980 10 15 0.0
> >> 16 1980 10 16 0.0
> >> 17 1980 10 17 41.2
> >> 18 1980 10 18 0.0
> >> 19 1980 10 19 0.0
> >> 20 1980 10 20 0.0
> >> 21 1980 10 21 0.0
> >> 22 1980 10 22 0.0
> >> 23 1980 10 23 0.0
> >> 24 1980 10 24 0.0
> >> 25 1980 10 25 0.0
> >> 26 1980 10 26 0.0
> >> 27 1980 10 27 2.0
> >> 28 1980 10 28 0.0
> >> 29 1980 10 29 0.0
> >> 30 1980 10 30 0.0
> >> 31 1980 10 31 0.0
> >> 32 1980 11 1 0.0
> >> 33 1980 11 2 0.0
> >> 34 1980 11 3 0.0
> >> 35 1980 11 4 0.0
> >> 36 1980 11 5 12.4
> >>
> >> the precipitation values extend from 1980 to 2005, but only for
october,
> >> november and december. I would like to plot just these 3 months for
the
> >> given time period (1980 - 2005). Is there a way to get these values in
> the
> >> x-axis. i.e. i need a plot with its axis reading "1980, 1981, 1982 ..
> 2005"
> >> (even if the months are specified then it should be of more use). for
> now,
> >> i get an axis like .. "0, 500, ... 2000" (the plot is giving the index
> >> values).
> >>
[[elided Yahoo spam]]
> >> ts.chn <- ts(chn.arr[1:2386,4], start=c(1980, 10), end=c(2005, 12),
> >> freq=365)
> >> plot(ts.chn)
> >>
> >> --
> >> Regards,
> >> Mahalakshmi
> >> Graduate Student
> >> #20, Department of Geography
> >> Michigan State University
> >> East Lansing, MI 48824
> >> [[alternative HTML version deleted]]
> >>
> >>
> >> _______________________________________________
> >> R-sig-Geo mailing list
> >> R-sig-Geo at stat.math.ethz.ch
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>
> >>
> >
>
[[alternative HTML version deleted]]
------------------------------
Message: 18
Date: Thu, 4 Nov 2010 16:05:11 -0700
From: Oliver Soong <osoong+r at gmail.com>
To: r-sig-geo at stat.math.ethz.ch
Subject: [R-sig-Geo] rgdal, multiband CCITT compression
Message-ID:
<AANLkTinVFAEBzfmjUZL84h7GMWpVjZtudfYrZ8dg264k at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
Hello,
I can't seem to write a multi-band 1-bit CCITT-compressed GeoTIFF
(using any of the CCITT compression methods). The file writes, but
when reloaded the data are wrong (I can't tell for sure whether this
is a read or write problem). I can write a 1-band 1-bit
CCITT-compressed GeoTIFF fine, and I can write a multi-band 1-bit
deflate-compressed GeoTIFF fine. This is on Windows XP x86, R 2.12.0,
rgdal 0.6-28, sp 0.9-72, GDAL 1.7.1, PROJ.4 4.7.1.
On a similar note, only GDAL seems to be able to handle any multi-band
1-bit GeoTIFFs it produces. ArcGIS 10 explicitly states (somewhere
deep in their documentation) that multi-band 1-bit TIFFs are not
supported. Even without compression, ENVI throws an error about an
array with a corrupted descriptor, then crashes out.
Is this a bug, or a limitation of TIFF, CCITT, or GDAL that I'm not aware
of?
Thoughts are appreciated,
Oliver
--
Oliver Soong
Donald Bren School of Environmental Science & Management
University of California, Santa Barbara
Santa Barbara, CA 93106-5131
805-893-7044 (office)
610-291-9706 (cell)
------------------------------
Message: 19
Date: Fri, 5 Nov 2010 10:18:34 +1100
From: Michael Sumner <mdsumner at gmail.com>
To: Oliver Soong <osoong+r at gmail.com>
Cc: r-sig-geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] rgdal, multiband CCITT compression
Message-ID:
<AANLkTi=SYUAuBQELQVesssJfrxXfqdLhdcwAt_L8kxNU at mail.gmail.com>
Content-Type: text/plain
Could you please provide an example of the problem with R code, so we have
something concrete to experiment with?
I could work out the required options and a sample dataset etc. but
presumably that is something you can do with less effort.
You should be able to investigate the generated file with the gdal command
line tools to figure out the problem.
Cheers, Mike.
On Fri, Nov 5, 2010 at 10:05 AM, Oliver Soong
<osoong+r at gmail.com<osoong%2Br at gmail.com>
> wrote:
> Hello,
>
> I can't seem to write a multi-band 1-bit CCITT-compressed GeoTIFF
> (using any of the CCITT compression methods). The file writes, but
> when reloaded the data are wrong (I can't tell for sure whether this
> is a read or write problem). I can write a 1-band 1-bit
> CCITT-compressed GeoTIFF fine, and I can write a multi-band 1-bit
> deflate-compressed GeoTIFF fine. This is on Windows XP x86, R 2.12.0,
> rgdal 0.6-28, sp 0.9-72, GDAL 1.7.1, PROJ.4 4.7.1.
>
> On a similar note, only GDAL seems to be able to handle any multi-band
> 1-bit GeoTIFFs it produces. ArcGIS 10 explicitly states (somewhere
> deep in their documentation) that multi-band 1-bit TIFFs are not
> supported. Even without compression, ENVI throws an error about an
> array with a corrupted descriptor, then crashes out.
>
> Is this a bug, or a limitation of TIFF, CCITT, or GDAL that I'm not aware
> of?
>
> Thoughts are appreciated,
> Oliver
>
>
> --
> Oliver Soong
> Donald Bren School of Environmental Science & Management
> University of California, Santa Barbara
> Santa Barbara, CA 93106-5131
> 805-893-7044 (office)
> 610-291-9706 (cell)
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Michael Sumner
Institute for Marine and Antarctic Studies, University of Tasmania
Hobart, Australia
e-mail: mdsumner at gmail.com
[[alternative HTML version deleted]]
------------------------------
Message: 20
Date: Thu, 4 Nov 2010 23:50:49 +0000
From: Craig Mundy <Craig.Mundy at utas.edu.au>
To: "r-sig-geo at stat.math.ethz.ch" <r-sig-geo at stat.math.ethz.ch>
Subject: [R-sig-Geo] Is there a function for counting polygons within
a multi-polygon SpatialPolygonDataFrame class
Message-ID:
<5504936889E8F947B474AFF792016B5D0353CC at MBXSBYN3.utas.ad.internal>
Content-Type: text/plain
Hi All,
I have searched for (but not found), a function for counting polygons within
each record of a multi-polygon SpatialPolygonDataFrame class? I had thought
of sifting through the slot information, but thought I would check if there
was a function I had overlooked before I embarked on that. I can easily do
this in SQL Server but would prefer to batch process in R along with some
other tasks.
Thanks
Craig
[[alternative HTML version deleted]]
------------------------------
Message: 21
Date: Thu, 4 Nov 2010 17:00:42 -0700
From: Oliver Soong <osoong+r at gmail.com>
To: Michael Sumner <mdsumner at gmail.com>
Cc: r-sig-geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] rgdal, multiband CCITT compression
Message-ID:
<AANLkTi=A33QZZMtn4vyC1ycH52MSd_rORKTU-LM5uR9+ at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
Here's an example:
require(rgdal)
formals(paste)$sep <- ""
SP27GTIF <- readGDAL(
system.file("pictures/SP27GTIF.TIF", package = "rgdal")[1],
output.dim=c(100,100)
)
img_1bit_1band <- SP27GTIF
img_1bit_1band at data$band1 <- as.integer(SP27GTIF at data$band1 > 85)
img_1bit_2band <- SP27GTIF
img_1bit_2band at data <- data.frame(
band1 = img_1bit_1band at data$band1,
band2 = as.integer(SP27GTIF at data$band1 > 170)
)
for(img in c("img_1bit_1band", "img_1bit_2band")) {
cat(paste("*** Correct Stats: ", img, "\n"))
print(sapply(get(img)@data, table))
for(compress in c("NONE", "DEFLATE", "LZW", "CCITTRLE",
"CCITTFAX3", "CCITTFAX4")) {
file <- paste(img, "_", compress, ".tif")
writeGDAL(get(img), file, type = "Byte",
options = c("NBITS=1", paste("COMPRESS=", compress)))
cat(paste("Compression: ", compress, " from ", file, "\n"))
reread <- readGDAL(file, silent = TRUE)
print(sapply(reread at data, table))
}
cat("\n")
}
The expected output is:
*** Correct Stats: img_1bit_1band
band1
0 3086
1 6914
Compression: NONE from img_1bit_1band_NONE.tif
band1
0 3086
1 6914
Compression: DEFLATE from img_1bit_1band_DEFLATE.tif
band1
0 3086
1 6914
Compression: LZW from img_1bit_1band_LZW.tif
band1
0 3086
1 6914
Compression: CCITTRLE from img_1bit_1band_CCITTRLE.tif
band1
0 3086
1 6914
Compression: CCITTFAX3 from img_1bit_1band_CCITTFAX3.tif
band1
0 3086
1 6914
Compression: CCITTFAX4 from img_1bit_1band_CCITTFAX4.tif
band1
0 3086
1 6914
*** Correct Stats: img_1bit_2band
band1 band2
0 3086 8353
1 6914 1647
Compression: NONE from img_1bit_2band_NONE.tif
band1 band2
0 3086 8353
1 6914 1647
Compression: DEFLATE from img_1bit_2band_DEFLATE.tif
band1 band2
0 3086 8353
1 6914 1647
Compression: LZW from img_1bit_2band_LZW.tif
band1 band2
0 3086 8353
1 6914 1647
Compression: CCITTRLE from img_1bit_2band_CCITTRLE.tif
band1 band2
0 6131 8783
1 3869 1217
Compression: CCITTFAX3 from img_1bit_2band_CCITTFAX3.tif
band1 band2
0 6131 8783
1 3869 1217
Compression: CCITTFAX4 from img_1bit_2band_CCITTFAX4.tif
band1 band2
0 6131 8783
1 3869 1217
Everything looks fine except CCITT compression with 2 bands.
I wasn't able to dig up Windows GDAL 1.7.1 binaries, and the 1.7.0
binaries from FW_Tools seem to give similar behavior, although it's a
little hard to check because neither ArcGIS nor ENVI seems to like any
of the multi-band files.
------------------------------
Message: 22
Date: Fri, 5 Nov 2010 11:14:51 +1100
From: Michael Sumner <mdsumner at gmail.com>
To: Oliver Soong <osoong+r at gmail.com>
Cc: r-sig-geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] rgdal, multiband CCITT compression
Message-ID:
<AANLkTin9vVU1T-_ucmPMxr2yOVpRWcxj5tfStXSHNT-j at mail.gmail.com>
Content-Type: text/plain
Ok, thanks. I can try later today with Manifold System, which may be. I'm
not sure of other non-GDAL tools to double check TIFFs with. (ArcGIS 10 uses
GDAL but I don't know to what extent).
Is the output above what you expect, or what you get? That's what I get.
Sorry to be dim, but I don't see where you've specified what's actually
wrong (and I haven't explored enough to see if that's easy to determine from
the code).
That's a cunning use of formals() to change defaults for a function - a bit
radical for my taste, but I like it.
Cheers, Mike.
On Fri, Nov 5, 2010 at 11:00 AM, Oliver Soong
<osoong+r at gmail.com<osoong%2Br at gmail.com>
> wrote:
> Here's an example:
>
>
>
>
> require(rgdal)
> formals(paste)$sep <- ""
>
> SP27GTIF <- readGDAL(
> system.file("pictures/SP27GTIF.TIF", package = "rgdal")[1],
> output.dim=c(100,100)
> )
> img_1bit_1band <- SP27GTIF
> img_1bit_1band at data$band1 <- as.integer(SP27GTIF at data$band1 > 85)
> img_1bit_2band <- SP27GTIF
> img_1bit_2band at data <- data.frame(
> band1 = img_1bit_1band at data$band1,
> band2 = as.integer(SP27GTIF at data$band1 > 170)
> )
>
> for(img in c("img_1bit_1band", "img_1bit_2band")) {
> cat(paste("*** Correct Stats: ", img, "\n"))
> print(sapply(get(img)@data, table))
> for(compress in c("NONE", "DEFLATE", "LZW", "CCITTRLE",
> "CCITTFAX3", "CCITTFAX4")) {
> file <- paste(img, "_", compress, ".tif")
> writeGDAL(get(img), file, type = "Byte",
> options = c("NBITS=1", paste("COMPRESS=", compress)))
> cat(paste("Compression: ", compress, " from ", file, "\n"))
> reread <- readGDAL(file, silent = TRUE)
> print(sapply(reread at data, table))
> }
> cat("\n")
> }
>
>
>
>
>
> The expected output is:
>
> *** Correct Stats: img_1bit_1band
> band1
> 0 3086
> 1 6914
> Compression: NONE from img_1bit_1band_NONE.tif
> band1
> 0 3086
> 1 6914
> Compression: DEFLATE from img_1bit_1band_DEFLATE.tif
> band1
> 0 3086
> 1 6914
> Compression: LZW from img_1bit_1band_LZW.tif
> band1
> 0 3086
> 1 6914
> Compression: CCITTRLE from img_1bit_1band_CCITTRLE.tif
> band1
> 0 3086
> 1 6914
> Compression: CCITTFAX3 from img_1bit_1band_CCITTFAX3.tif
> band1
> 0 3086
> 1 6914
> Compression: CCITTFAX4 from img_1bit_1band_CCITTFAX4.tif
> band1
> 0 3086
> 1 6914
>
> *** Correct Stats: img_1bit_2band
> band1 band2
> 0 3086 8353
> 1 6914 1647
> Compression: NONE from img_1bit_2band_NONE.tif
> band1 band2
> 0 3086 8353
> 1 6914 1647
> Compression: DEFLATE from img_1bit_2band_DEFLATE.tif
> band1 band2
> 0 3086 8353
> 1 6914 1647
> Compression: LZW from img_1bit_2band_LZW.tif
> band1 band2
> 0 3086 8353
> 1 6914 1647
> Compression: CCITTRLE from img_1bit_2band_CCITTRLE.tif
> band1 band2
> 0 6131 8783
> 1 3869 1217
> Compression: CCITTFAX3 from img_1bit_2band_CCITTFAX3.tif
> band1 band2
> 0 6131 8783
> 1 3869 1217
> Compression: CCITTFAX4 from img_1bit_2band_CCITTFAX4.tif
> band1 band2
> 0 6131 8783
> 1 3869 1217
>
>
>
> Everything looks fine except CCITT compression with 2 bands.
>
> I wasn't able to dig up Windows GDAL 1.7.1 binaries, and the 1.7.0
> binaries from FW_Tools seem to give similar behavior, although it's a
> little hard to check because neither ArcGIS nor ENVI seems to like any
> of the multi-band files.
>
--
Michael Sumner
Institute for Marine and Antarctic Studies, University of Tasmania
Hobart, Australia
e-mail: mdsumner at gmail.com
[[alternative HTML version deleted]]
------------------------------
Message: 23
Date: Fri, 5 Nov 2010 11:15:27 +1100
From: Michael Sumner <mdsumner at gmail.com>
To: Oliver Soong <osoong+r at gmail.com>
Cc: r-sig-geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] rgdal, multiband CCITT compression
Message-ID:
<AANLkTi=7G1dK-o-_yqOOc5KHyNjJRq-c5NgYFe9hMAWA at mail.gmail.com>
Content-Type: text/plain
Sorry, I meant:
"I can try later today with Manifold System, which may be helpful as an
independent tool, maybe".
On Fri, Nov 5, 2010 at 11:14 AM, Michael Sumner <mdsumner at gmail.com> wrote:
> Ok, thanks. I can try later today with Manifold System, which may be. I'm
> not sure of other non-GDAL tools to double check TIFFs with. (ArcGIS 10
uses
> GDAL but I don't know to what extent).
>
> Is the output above what you expect, or what you get? That's what I get.
> Sorry to be dim, but I don't see where you've specified what's actually
> wrong (and I haven't explored enough to see if that's easy to determine
from
> the code).
>
> That's a cunning use of formals() to change defaults for a function - a
bit
> radical for my taste, but I like it.
>
> Cheers, Mike.
>
> On Fri, Nov 5, 2010 at 11:00 AM, Oliver Soong
<osoong+r at gmail.com<osoong%2Br at gmail.com>
> > wrote:
>
>> Here's an example:
>>
>>
>>
>>
>> require(rgdal)
>> formals(paste)$sep <- ""
>>
>> SP27GTIF <- readGDAL(
>> system.file("pictures/SP27GTIF.TIF", package = "rgdal")[1],
>> output.dim=c(100,100)
>> )
>> img_1bit_1band <- SP27GTIF
>> img_1bit_1band at data$band1 <- as.integer(SP27GTIF at data$band1 > 85)
>> img_1bit_2band <- SP27GTIF
>> img_1bit_2band at data <- data.frame(
>> band1 = img_1bit_1band at data$band1,
>> band2 = as.integer(SP27GTIF at data$band1 > 170)
>> )
>>
>> for(img in c("img_1bit_1band", "img_1bit_2band")) {
>> cat(paste("*** Correct Stats: ", img, "\n"))
>> print(sapply(get(img)@data, table))
>> for(compress in c("NONE", "DEFLATE", "LZW", "CCITTRLE",
>> "CCITTFAX3", "CCITTFAX4")) {
>> file <- paste(img, "_", compress, ".tif")
>> writeGDAL(get(img), file, type = "Byte",
>> options = c("NBITS=1", paste("COMPRESS=", compress)))
>> cat(paste("Compression: ", compress, " from ", file, "\n"))
>> reread <- readGDAL(file, silent = TRUE)
>> print(sapply(reread at data, table))
>> }
>> cat("\n")
>> }
>>
>>
>>
>>
>>
>> The expected output is:
>>
>> *** Correct Stats: img_1bit_1band
>> band1
>> 0 3086
>> 1 6914
>> Compression: NONE from img_1bit_1band_NONE.tif
>> band1
>> 0 3086
>> 1 6914
>> Compression: DEFLATE from img_1bit_1band_DEFLATE.tif
>> band1
>> 0 3086
>> 1 6914
>> Compression: LZW from img_1bit_1band_LZW.tif
>> band1
>> 0 3086
>> 1 6914
>> Compression: CCITTRLE from img_1bit_1band_CCITTRLE.tif
>> band1
>> 0 3086
>> 1 6914
>> Compression: CCITTFAX3 from img_1bit_1band_CCITTFAX3.tif
>> band1
>> 0 3086
>> 1 6914
>> Compression: CCITTFAX4 from img_1bit_1band_CCITTFAX4.tif
>> band1
>> 0 3086
>> 1 6914
>>
>> *** Correct Stats: img_1bit_2band
>> band1 band2
>> 0 3086 8353
>> 1 6914 1647
>> Compression: NONE from img_1bit_2band_NONE.tif
>> band1 band2
>> 0 3086 8353
>> 1 6914 1647
>> Compression: DEFLATE from img_1bit_2band_DEFLATE.tif
>> band1 band2
>> 0 3086 8353
>> 1 6914 1647
>> Compression: LZW from img_1bit_2band_LZW.tif
>> band1 band2
>> 0 3086 8353
>> 1 6914 1647
>> Compression: CCITTRLE from img_1bit_2band_CCITTRLE.tif
>> band1 band2
>> 0 6131 8783
>> 1 3869 1217
>> Compression: CCITTFAX3 from img_1bit_2band_CCITTFAX3.tif
>> band1 band2
>> 0 6131 8783
>> 1 3869 1217
>> Compression: CCITTFAX4 from img_1bit_2band_CCITTFAX4.tif
>> band1 band2
>> 0 6131 8783
>> 1 3869 1217
>>
>>
>>
>> Everything looks fine except CCITT compression with 2 bands.
>>
>> I wasn't able to dig up Windows GDAL 1.7.1 binaries, and the 1.7.0
>> binaries from FW_Tools seem to give similar behavior, although it's a
>> little hard to check because neither ArcGIS nor ENVI seems to like any
>> of the multi-band files.
>>
>
>
>
> --
> Michael Sumner
> Institute for Marine and Antarctic Studies, University of Tasmania
> Hobart, Australia
> e-mail: mdsumner at gmail.com
>
--
Michael Sumner
Institute for Marine and Antarctic Studies, University of Tasmania
Hobart, Australia
e-mail: mdsumner at gmail.com
[[alternative HTML version deleted]]
------------------------------
Message: 24
Date: Fri, 5 Nov 2010 11:33:43 +1100
From: Michael Sumner <mdsumner at gmail.com>
To: Craig Mundy <Craig.Mundy at utas.edu.au>
Cc: "r-sig-geo at stat.math.ethz.ch" <r-sig-geo at stat.math.ethz.ch>
Subject: Re: [R-sig-Geo] Is there a function for counting polygons
within a multi-polygon SpatialPolygonDataFrame class
Message-ID:
<AANLkTimgrHqFNSwqns3_Y_BFYR_LHyrX0hf5KYq-jjgF at mail.gmail.com>
Content-Type: text/plain; charset=UTF-8
This will give you the number of separate "rings", but that won't
allow you to differentiate between separate simple polygons, and rings
that represent "holes":
sapply(SrDf at polygons, function(x) length(slot(x, "Polygons")))
It needs more work to check each individual "Polygon" in "x" to see if
slot(x[[i]], "hole") is TRUE. I think this will do it:
sapply(SrDf at polygons, function(x) sum( !sapply(slot(x, "Polygons"),
slot, "hole")))
Note that use slot(x, "name") is more verbose that x at name, but either
is not recommended for normal use and should be bundled into methods,
or helper functions.
I think it's worth delving into the structure of these objects, as
it's fairly simple once you get used to it and writing summary
functions is not hard.
Cheers, Mike.
On Fri, Nov 5, 2010 at 10:50 AM, Craig Mundy <Craig.Mundy at utas.edu.au>
wrote:
>
> Hi All,
>
> I have searched for (but not found), a function for counting polygons
within each record of a multi-polygon SpatialPolygonDataFrame class? I had
thought of sifting through the slot information, but thought I would check
if there was a function I had overlooked before I embarked on that. I can
easily do this in SQL Server but would prefer to batch process in R along
with some other tasks.
>
> Thanks
>
> Craig
>
> ? ? ? ?[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Michael Sumner
Institute for Marine and Antarctic Studies, University of Tasmania
Hobart, Australia
e-mail: mdsumner at gmail.com
------------------------------
Message: 25
Date: Thu, 4 Nov 2010 17:40:27 -0700
From: Oliver Soong <osoong+r at gmail.com>
To: Michael Sumner <mdsumner at gmail.com>
Cc: r-sig-geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] rgdal, multiband CCITT compression
Message-ID:
<AANLkTimcNKFC4exW1KjYyJq1DxR1aTka+VqUbTvquosk at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
Sorry, the word "expected" was a bad choice. The output is what I
get. I'm highlighting two tables that should be the same, but aren't.
I appreciate you taking the time to poke around with me.
*** Correct Stats: img_1bit_2band
band1 band2
0 3086 8353
1 6914 1647
Compression: CCITTRLE from img_1bit_2band_CCITTRLE.tif
band1 band2
0 6131 8783
1 3869 1217
On Thu, Nov 4, 2010 at 5:14 PM, Michael Sumner <mdsumner at gmail.com> wrote:
> Ok, thanks. I can try later today with Manifold System, which may be. I'm
> not sure of other non-GDAL tools to double check TIFFs with. (ArcGIS 10
uses
> GDAL but I don't know to what extent).
>
> Is the output above what you expect, or what you get? That's what I get.
> Sorry to be dim, but I don't see where you've specified what's actually
> wrong (and I haven't explored enough to see if that's easy to determine
from
> the code).
>
> That's a?cunning use of formals() to?change defaults for a function?- a
bit
> radical for my taste, but I like it.
>
> Cheers, Mike.
>
> On Fri, Nov 5, 2010 at 11:00 AM, Oliver Soong <osoong+r at gmail.com> wrote:
>>
>> Here's an example:
>>
>>
>>
>>
>> require(rgdal)
>> formals(paste)$sep <- ""
>>
>> SP27GTIF <- readGDAL(
>> ?system.file("pictures/SP27GTIF.TIF", package = "rgdal")[1],
>> ?output.dim=c(100,100)
>> )
>> img_1bit_1band <- SP27GTIF
>> img_1bit_1band at data$band1 <- as.integer(SP27GTIF at data$band1 > 85)
>> img_1bit_2band <- SP27GTIF
>> img_1bit_2band at data <- data.frame(
>> ?band1 = img_1bit_1band at data$band1,
>> ?band2 = as.integer(SP27GTIF at data$band1 > 170)
>> )
>>
>> for(img in c("img_1bit_1band", "img_1bit_2band")) {
>> ?cat(paste("*** Correct Stats: ", img, "\n"))
>> ?print(sapply(get(img)@data, table))
>> ?for(compress in c("NONE", "DEFLATE", "LZW", "CCITTRLE",
>> ? ? ? ? ? ? ? ? ? ?"CCITTFAX3", "CCITTFAX4")) {
>> ? ?file <- paste(img, "_", compress, ".tif")
>> ? ?writeGDAL(get(img), file, type = "Byte",
>> ? ? ?options = c("NBITS=1", paste("COMPRESS=", compress)))
>> ? ?cat(paste("Compression: ", compress, " from ", file, "\n"))
>> ? ?reread <- readGDAL(file, silent = TRUE)
>> ? ?print(sapply(reread at data, table))
>> ?}
>> ?cat("\n")
>> }
>>
>>
>>
>>
>>
>> The expected output is:
>>
>> *** Correct Stats: img_1bit_1band
>> ?band1
>> 0 ?3086
>> 1 ?6914
>> Compression: NONE from img_1bit_1band_NONE.tif
>> ?band1
>> 0 ?3086
>> 1 ?6914
>> Compression: DEFLATE from img_1bit_1band_DEFLATE.tif
>> ?band1
>> 0 ?3086
>> 1 ?6914
>> Compression: LZW from img_1bit_1band_LZW.tif
>> ?band1
>> 0 ?3086
>> 1 ?6914
>> Compression: CCITTRLE from img_1bit_1band_CCITTRLE.tif
>> ?band1
>> 0 ?3086
>> 1 ?6914
>> Compression: CCITTFAX3 from img_1bit_1band_CCITTFAX3.tif
>> ?band1
>> 0 ?3086
>> 1 ?6914
>> Compression: CCITTFAX4 from img_1bit_1band_CCITTFAX4.tif
>> ?band1
>> 0 ?3086
>> 1 ?6914
>>
>> *** Correct Stats: img_1bit_2band
>> ?band1 band2
>> 0 ?3086 ?8353
>> 1 ?6914 ?1647
>> Compression: NONE from img_1bit_2band_NONE.tif
>> ?band1 band2
>> 0 ?3086 ?8353
>> 1 ?6914 ?1647
>> Compression: DEFLATE from img_1bit_2band_DEFLATE.tif
>> ?band1 band2
>> 0 ?3086 ?8353
>> 1 ?6914 ?1647
>> Compression: LZW from img_1bit_2band_LZW.tif
>> ?band1 band2
>> 0 ?3086 ?8353
>> 1 ?6914 ?1647
>> Compression: CCITTRLE from img_1bit_2band_CCITTRLE.tif
>> ?band1 band2
>> 0 ?6131 ?8783
>> 1 ?3869 ?1217
>> Compression: CCITTFAX3 from img_1bit_2band_CCITTFAX3.tif
>> ?band1 band2
>> 0 ?6131 ?8783
>> 1 ?3869 ?1217
>> Compression: CCITTFAX4 from img_1bit_2band_CCITTFAX4.tif
>> ?band1 band2
>> 0 ?6131 ?8783
>> 1 ?3869 ?1217
>>
>>
>>
>> Everything looks fine except CCITT compression with 2 bands.
>>
>> I wasn't able to dig up Windows GDAL 1.7.1 binaries, and the 1.7.0
>> binaries from FW_Tools seem to give similar behavior, although it's a
>> little hard to check because neither ArcGIS nor ENVI seems to like any
>> of the multi-band files.
>
>
>
> --
> Michael Sumner
> Institute for Marine and Antarctic Studies, University of Tasmania
> Hobart, Australia
> e-mail: mdsumner at gmail.com
>
--
Oliver Soong
Donald Bren School of Environmental Science & Management
University of California, Santa Barbara
Santa Barbara, CA 93106-5131
805-893-7044 (office)
610-291-9706 (cell)
------------------------------
Message: 26
Date: Fri, 5 Nov 2010 07:15:32 -0200
From: "Rodrigo Aluizio" <r.aluizio at gmail.com>
To: "R Help" <r-sig-geo at stat.math.ethz.ch>
Subject: [R-sig-Geo] Defining Holes Inside a Polygon
Message-ID: <000001cb7cca$02076c30$06164490$@gmail.com>
Content-Type: text/plain; charset="iso-8859-1"
Hi List.
I?m trying to solve an issue here, but it?s getting tricky.
I have an MapGen .dat file (attached) that contains about 39 closed lines.
Three of them (the largest ones) will be my main polygons and the rest of
them must become holes inside one of the three major areas. I?m able to
transform all the closed lines into SpatialPolygons and even create
shapefiles with them, but how can I modify the ?hole? slot of the smaller
polygons so I can keep them transparent when overlaying another map (those
holes will overlay islands of a raster image)?
Thank you for the attention and help.
Regards.
-------------------------------------------------------------------
MSc. <mailto:r.aluizio at gmail.com> Rodrigo Aluizio
Centro de Estudos do Mar/UFPR
Laborat?rio de Foramin?feros e Micropaleontologia Ambiental
Avenida Beira Mar s/n - CEP 83255-971
Pontal do Paran? - PR - Brazil
-------------- next part --------------
An HTML attachment was scrubbed...
URL:
<https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20101105/1d7fee34/atta
chment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: PCForams.dat
Type: application/octet-stream
Size: 133962 bytes
Desc: not available
URL:
<https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20101105/1d7fee34/atta
chment-0001.obj>
------------------------------
Message: 27
Date: Fri, 5 Nov 2010 11:01:05 +0100 (CET)
From: Roger Bivand <Roger.Bivand at nhh.no>
To: Rodrigo Aluizio <r.aluizio at gmail.com>
Cc: R Help <r-sig-geo at stat.math.ethz.ch>
Subject: Re: [R-sig-Geo] Defining Holes Inside a Polygon
Message-ID: <alpine.LRH.2.00.1011051048370.27156 at reclus.nhh.no>
Content-Type: text/plain; charset="iso-8859-1"; Format="flowed"
On Fri, 5 Nov 2010, Rodrigo Aluizio wrote:
> Hi List.
>
> I?m trying to solve an issue here, but it?s getting tricky.
>
> I have an MapGen .dat file (attached) that contains about 39 closed lines.
> Three of them (the largest ones) will be my main polygons and the rest of
> them must become holes inside one of the three major areas. I?m able to
> transform all the closed lines into SpatialPolygons and even create
> shapefiles with them, but how can I modify the ?hole? slot of the smaller
> polygons so I can keep them transparent when overlaying another map (those
> holes will overlay islands of a raster image)?
>
Next time avoid sending 200K to over 1700 people, please, we have a planet
to protect! Put any desired attachments on a website and provide a link.
You should also have provided the code (in your message) to get from:
library(maptools)
l1 <- MapGen2SL("PCForams.dat")
to your problem - as those choices may be the cause of the problem; you
have not provided the output of sessionInfo() either. Once things are
lined up, you should - if you may - say:
gpclibPermit() # very restrictive license
pls_out <- lapply(slot(SpPols0, "polygons"), checkPolygonsHoles)
SpPols1 <- SpatialPolygons(pls_out)
Or if you can install rgeos from source on R-Forge, equivalently:
library(rgeos)
SpPols1 <- createSPComment(SpPols0)
Hope this helps,
Roger
>
>
> Thank you for the attention and help.
>
> Regards.
>
>
>
> -------------------------------------------------------------------
>
> MSc. <mailto:r.aluizio at gmail.com> Rodrigo Aluizio
>
> Centro de Estudos do Mar/UFPR
> Laborat?rio de Foramin?feros e Micropaleontologia Ambiental
> Avenida Beira Mar s/n - CEP 83255-971
> Pontal do Paran? - PR - Brazil
>
>
>
>
--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no
------------------------------
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
End of R-sig-Geo Digest, Vol 87, Issue 5
More information about the R-sig-Geo
mailing list