On Wed, 10 Feb 2010, Agustin Lobo wrote:
FYI This problem was actually identified because I was getting wrong results at calculating averaged raster values for intersecting line segments. I include the last message sent by Robert. His trick is a good workaround by now.
And there is a similar mechanism in sp. However, please find out and report to the plugin author why the plugin is flipping GTiffs.
Roger
Agus ---------- Forwarded message ---------- From: "Robert J. Hijmans" <r.hijm...@gmail.com> Date: Wed, 10 Feb 2010 09:01:30 -0800 Subject: Re: [R-sig-Geo] Overlay spatial lines on raster To: Agustin Lobo <agustin.l...@ija.csic.es> Ah! The fundamentals seems to be beyond my reach (in (r)gdal) but the simple raster patch if you *know* this is happening is: r = flip(r) Robert On Wed, Feb 10, 2010 at 4:19 AM, Agustin Lobo <agustin.l...@ija.csic.es> wrote:Robert, Sorry I could not react earlier to this mystery, I've found a clue now: the raster is reversed in R in respect to QGIS (see attachment). I've uploaded test2.tif (identical to test.tif but with the appropriate projection info) to http://sites.google.com/site/eospansite/dummy/test2.tif The same problem arises with raster() and with readGDAL(), so we should CC to sig-geo at some point? This is what I do:r2 <- raster("test2.tif") plot(r2) GDALinfo("test2.tif")rows 256 columns 256 bands 1 origin.x 424389 origin.y 4636016 res.x 345.0078 res.y 194.0195 oblique.x 0 oblique.y 0 driver GTiff projection +proj=utm +zone=31 +ellps=intl +units=m +no_defs file test2.tif apparent band summary: GDType Bmin Bmax 1 Float64 -4294967295 4294967295 Metadata: AREA_OR_POINT=Arear3 <- readGDAL("test2.tif")test2.tif has GDAL driver GTiff and has 256 rows and 256 columnsimage(r3)I'm not sure I'm not doing something really stupid here, I apologize in that case! Thanks Agus 2010/2/1, Robert J. Hijmans <r.hijm...@gmail.com>:Agus, I am not sure what is going on here. But this example seems to work for me: library(raster) setwd('c:/downloads') load("xter.rda") r <- raster("test.tif") projection(r) <- projection(xter) # to speed things up for testing r = crop(r, extent(423517, 454000, 4675610 , 4686260 )) lr = linesToRaster(xter, r, field = "ID_GRAFIC", progress='text') z = zonal(r, lr, mean) z$mean = as.numeric(z$mean) # this will not be necessary in the next version d <- x...@data dd <- merge(d, z, by.x="ID_GRAFIC", by.y="zone", sort=F, all.x=T) x...@data <- dd[order(dd[,1]),] # just for testing test = linesToRaster(xter, r, field = "mean", progress='text') par(mfrow=c(1,3)) plot(lr) plot(r) lines(xter) plot(test) # looks good, also in spplot x11() spplot(xter, 'mean') writeOGR(xter, dsn="xter",layer="xtertest",driver="ESRI Shapefile") On Mon, Feb 1, 2010 at 1:21 AM, Agustin Lobo <agustin.l...@ija.csic.es> wrote:Robert, Thanks. The mechanics work, but I get inconsistent results. This is what I do (there is also a minor problem with writeOGR() that I'm circumventing by now): lr = linesToRaster(xter, lr, field = "ID_GRAFIC" ) xterpdata2 = zonal(r, lr, mean) xterpdata2 = xterpdata2[order(xterpdata2[,1]),] xterpdata3 <- na.omit(xterpdata2) dim(xterpdata3) dim(xterpdata2) xter2 <- xter xter2data <- xt...@data xter2data <- merge(xter2data,xterpdata3,by.x="ID_GRAFIC",by.y="zone",sort=F,all.x=T) xter2data <- xter2data[order(xter2data[,1]),] xt...@data <- xter2data writeOGR(xter2, dsn="xter2",layer="xter2",driver="ESRI Shapefile") #and get an error: #Error in writeOGR(xter2, dsn = "xter2", layer = "xter2", driver = "ESRI Shapefile") : # Can't convert columns of class: array; column names: mean #class(xt...@data$mean) #[1] "array" #which I curcunvent by: xt...@data$mean <- as.numeric(xt...@data$mean) writeOGR(xter2, dsn="xter2",layer="xter2",driver="ESRI Shapefile") writeRaster(lr,filename="lr",format="GTiff") Now, we check for example the result of zonal()xterpdata2[xterpdata2[,1]==8500,]zone mean 8500 8500 4.253456 But in the QGIS we see that the values of test.tif (r in R) are very different (~19 in test (see to the left), while field "mean" is 4.253, see attached picture). Could be the problem is in zonal() ? Agus Robert J. Hijmans wrote:Hola Agus, Thanks for sending these data. That was very helpful. This is a whole different (subsetting) problem for line segments consisting of two points only (with simple solution, "drop=FALSE"). Nice to have found it. Your example now runs for me, although it takes three minutes on my PC (I'll speed it up eventually, let's first make sure it works). You can use the attached updated package (or wait till tomorrow for R-Forge). Robert On Sun, Jan 31, 2010 at 3:18 PM, Agustin Lobo <agustin.l...@ija.csic.es> wrote:Thanks, actually the sourceforge version was ready when I tried again (few minutes ago) Still problems. I've uploaded the files for you to try if you wish: http://sites.google.com/site/eospansite/dummy/xter.rda http://sites.google.com/site/eospansite/dummy/test.tif It would be great if you could run: load("xter.rda") r <- raster("test.tif") projection(r) <- projection(xter) lr = linesToRaster(xter, r, field = "ID_GRAFIC" ) I get: Error in subset.matrix(xyxy, !((xyxy[, 2] > maxy & xyxy[, 4] > maxy) | : subscript out of bounds Calls: linesToRaster -> .getCols -> subset -> subset.matrix Agus Robert J. Hijmans wrote:Hi Agus, Attached is the linux (is that your o.s.?) package for local install. Robert On Sun, Jan 31, 2010 at 2:00 AM, Agustin Lobo <alobolis...@gmail.com> wrote:Looks great, is there any way I can try right now? Otherwise I'll wait until 0.9.9-5 is ready for installation. I tried now and only 0.9.9-4 is available. Thanks! Agus Robert J. Hijmans wrote:Hi Agus, thank you. You were right about the cause of the problem --- I patched it in version 0.9.9-5, now on R-Forge (and should be available for automatic install within 24 hours); the function could use some optimization for speed but at least it works now for SpatialLines* with a larger extent then the target RasterLayer. # simple example library(raster) library(maptools) data(wrld_simpl) r = raster(xmn=-10, xmx=20, ymn=37, ymx=43, ncol=360, nrow=72) r = linesToRaster(wrld_simpl, r, progress='text') plot(r) plot(wrld_simpl, add=T) # example using zonal as I proposed in your original question, but here with polygonsToRaster # What is (roughly) the mean latitude of each country? x = raster() # res(x) = 0.1 # higher res, slower, but more small countries are included x = polygonsToRaster(wrld_simpl, x, progress='text') y = raster(x) y[] = rep(yFromRow(y, 1:nrow(y)), each=ncol(y)) z = zonal(y, x, mean) res= data.frame(wrld_simpl[z[,'zone'], 'NAME'], lat=z[,2]) res[order(res[,2]), ] Robert On Sat, Jan 30, 2010 at 5:03 AM, Agustin Lobo <alobolis...@gmail.com> wrote:Thanks Robert, The raster package is a magician's chest! I get an error, though:class(xter)[1] "SpatialLinesDataFrame" attr(,"package") [1] "sp"class(r)[1] "RasterLayer" attr(,"package") [1] "raster"lr = linesToRaster(xter, r, field = "ID_GRAFIC" )[1] "A" rows cols [1,] 1 -1 [2,] -1 -1 [,1] [,2] [1,] 413809.3 4685685 [1] 413786.5 4685646.9 413810.1 4685686.4 Error in col1:col2 : NA/NaN argument Calls: linesToRaster -> .getCols I think because the extent of the raster is smaller than the extent of the linesextent(xter)class : Extent xmin : 412792 xmax : 517648.1 ymin : 4624794 ymax : 4698657extent(r)class : Extent xmin : 424389 xmax : 512711 ymin : 4636016 ymax : 4685685 which I can visualize with:plot(extent(xter)) plot(extent(r),col="red",add=T)The problem is that I cannot find the way to clip xter to the extent of r. pruneMap() assumes a map object and I guess that we cannot convert from Spatial Lines to map. I've trieddelme <- xter de...@bbox <- bbox(r)but still get the same problem:lr = linesToRaster(delme, r, field = "ID_GRAFIC" )[1] "A" rows cols [1,] 1 -1 [2,] -1 -1 [,1] [,2] [1,] NA NA [1] 413786.5 4685646.9 413810.1 4685686.4 Error in col1:col2 : NA/NaN argument Calls: linesToRaster -> .getCols Agus 2010/1/29 Robert J. Hijmans <r.hijm...@gmail.com>How about something like: library(raster) r = linesToRaster(xter, xterpdata, field = "ID" ) v = zonal(xterpdata, r, mean) v now shold has the mean value of the raster for each line ID and you can merge those back On Fri, Jan 29, 2010 at 9:54 AM, Agustin Lobo <alobolis...@gmail.com> wrote:Partially answering myself: This is what I've done until now: xter <- readOGR(dsn="/media/Transcend/PROVI/MARICEL2/x_terRiverBasin",layer="x_ter")class(xter)[1] "SpatialLinesDataFrame" attr(,"package") [1] "sp" xterp <- coordinates(xter) xterp <- sapply(xterp, function(x) do.call("rbind", x)) xterp <- do.call("rbind",xterp) as told by Roger. Now: require(raster) r <- raster("test.tif") xterpdata <- xyValues(r,xterp) xterp2 <- data.frame(utmx=xterp[,1], utmy=xterp[,2], NO3_2007=xterpdata) coordinates(xterp2) <- c("utmx", "utmy") So I get to step 2, but don't see how to convert from the points to the lines again (to get xterl). The problem is that when I did: xterp <- sapply(xterp, function(x) do.call("rbind", x)) xterp <- do.call("rbind",xterp) I lost the IDs of the lines. So I would need to keep the IDs of the lines for each point in the result of coordinates() to be able of going back to the lines once I've put the raster values in the data slot of the xterp I'll keep on, but any help would be much appreciated. Agus 2010/1/29 Agustin Lobo <alobolis...@gmail.com>Hi! My goal is to colorcode a river network according to the values of a raster map (which comes from an interpolation). The spatial lines are imported from a shapefile. The raster could be imported or I could run the interpolation in R. As what I'm planing could be rather involved, I prefer asking before: 1. Convert spatial lines to spatial points as discussed earlier this week. 2. Overlay points on the raster and get the values for the points in a Spatial Points Data Frame 3. Convert the Spatial Points Data Frame into Spatial Lines Data Frame 4. writeOGR() to shapefile. The main question is: is the data frame of the points going to be transferred to the lines at step 3? THnaks Agus[[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo-- Dr. Agustin Lobo Institut de Ciencies de la Terra "Jaume Almera" (CSIC) LLuis Sole Sabaris s/n 08028 Barcelona Spain Tel. 34 934095410 Fax. 34 934110012 email: agustin.l...@ija.csic.es http://www.ija.csic.es/gt/obster-- Dr. Agustin Lobo Institut de Ciencies de la Terra "Jaume Almera" (CSIC) LLuis Sole Sabaris s/n 08028 Barcelona Spain Tel. 34 934095410 Fax. 34 934110012 email: agustin.l...@ija.csic.es http://www.ija.csic.es/gt/obster-- Dr. Agustin Lobo Institut de Ciencies de la Terra "Jaume Almera" (CSIC) Lluis Sole Sabaris s/n 08028 Barcelona Spain Tel. 34 934095410 Fax. 34 934110012 e-mail agustin.l...@ija.csic.es http://www.ija.csic.es/gt/obster-- Dr. Agustin Lobo Institut de Ciencies de la Terra "Jaume Almera" (CSIC) Lluis Sole Sabaris s/n 08028 Barcelona Spain Tel. 34 934095410 Fax. 34 934110012 e-mail agustin.l...@ija.csic.es http://www.ija.csic.es/gt/obster _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
-- 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.biv...@nhh.no
_______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
