Lorenzo Alfieri | 22 Jul 16:04 2014
Picon

apply interpolation functions on a map

Dear list,
I can't find a solution to this problem:
Let's say I have a map of water depth due to a flood , and the corresponding land-use:

library(raster)
WaterDepth<-raster(matrix(1:20, nrow = 5, ncol = 4))
LandUse<-raster(WaterDepth)
LandUse[]<-c(rep(1,12),rep(2,8))

# each land use class has a different relation between water depth and damage, which is defined per points
and that I can interpolate linearly:
Damage<-data.frame(V1=c(0,100,150,230,340,340),V2=c(0,55,120,160,225,225))
Depths<-c(0,0.5,1,2,4,100)     #Water depth in meters
DDFun<-lapply(Damage, approxfun, x=Depths)

# Now if I want to know the estimated cost of the damage on a specific grid point it seems to work:
i<-3
Cost<-raster(WaterDepth)
Cost[i] <- DDFun[[LandUse[i]]](WaterDepth[i])

> Cost[i]

340 

But it doesn't work if I try to apply it to the entire map extent:

> Cost[] <- DDFun[[LandUse[]]](WaterDepth[])
Error in DDFun[[LandUse[]]] : recursive indexing failed at level 2

Any suggestion on how to do this (possibly without looping on each land use class)?
(Continue reading)

Dominik Schneider | 22 Jul 02:26 2014
Picon

extracting raster data with point and polygons

Hi -
I have a 2 meter DEM, a 30m satellite image and some point data. Basically,
I would like to extract the 2m dem pixels that correspond with the 30m
satellite data, but only for the satellite pixels that correspond with a
point observation. i've been trying to get at the answer from the sp
vignette 'over' but am stuck.
an example:
#define DEM
dem2m=raster(x=matrix(runif(22500,min=2000,max=2500),nrow=150),crs='+proj=utm
+zone=13 +datum=NAD83')
extent(dem2m)=c(0,300,0,3000)
#
#define satellite image
fsca.dem=raster(x=matrix(runiff(100,0,100),nrow=10),crs='+proj=utm +zone=13
+datum=NAD83')
extent(dem2m)=c(0,300,0,3000)
#
#define point data
obs=data.frame(x=runif(5,0,300),y=runif(5,0,300),sd=runif(5,0,5))
coordinates(obs)=~x+y
proj4string(obs)='+proj=utm +zone=13 +datum=NAD83'
#
#--- in my case, I resampled the sat image to the dem grid. but thats
already taken care of in the example.
##dem.pe=projectExtent(dem2m,crs=projection(dem2m))
##res(dem.pe) <- 30
##fsca.dem=projectRaster(fsca,dem.pe,method='ngb')
#
# This is where I am lost
# ----------- I thought I'd use polygons to extract the dem data. The catch
(Continue reading)

Roland Pape | 21 Jul 21:22 2014
Picon

Best way to provide meaningful interpretation of axes out of GNESFA, KSelect and Niche analyses

Dear list,

while working with the adehabitatHS package to explore habitat use of different animal species I wondered
about the best solution to provide a meaningful interpretation of the axes resulting from the various
approaches. 
According to the vignette of the package (p 19), it is preferred to give meaning to the axes of the Gnesfa with
the help of correlations between the environmental variables and the axes of the Gnesfa rather than using
the coefficients of the variables. The reason to do so (correlated variables may render the
interpretation of coefficients difficult) holds also for other approaches like the Kselect and Niche
analyses. However, their output provides per default only variable coordinates and normed scores with
the latter being used to give meaning to the axes. Wouldn't it be appropriate to use correlations here as
well and if yes, how could they be achieved?

Kind regards,
Roland
Edoardo Baldoni | 21 Jul 12:27 2014
Picon

Logit spatial predictions outside range [0,1]

Hello R-sig-Geo,

my name is Edoardo and this is the first time I am asking for help here.

I am going through the "sdm" vignette of package "dismo" and, when using
the Logistic regression model of presence/absence(background) data of
bradypus (p. 53) for spatial prediction, I get values outside the range
[0,1]. I think this should not happen as using the Logit cumulative
distribution function should limit them in the [0,1] range.

It follows the code with the data:

## This is the dataset I use. It is just a random sample of the original
envtrain dataset of
## the "sdm" vignette

> envtrain0
    pa bio1 bio12 bio16 bio17 bio5 bio6 bio7 bio8 biome
735  0  265  2618   954   391  314  224   90  264     1
634  0  264  2590  1007   203  327  200  127  263     1
861  0  243  1783   916    29  335  132  202  247     1
70   1  261  2791  1077   148  310  218   92  254     1
48   1  257  3575  1295   359  316  202  114  258     1
511  0  220  1583   457   281  339  114  225  226     7
344  0  230  3352  1436   247  294  167  127  228     1
453  0  242   240   140    14  308  186  121  260    13
325  0  240  1541   447   302  291  186  104  249     1
193  0  232  1130   506    87  311  123  188  256     1
568  0  257  2269   774   246  317  188  128  259     1
664  0  190  1512   544   251  282   85  196  219     1
(Continue reading)

Kátia Emidio | 20 Jul 17:51 2014
Picon

Local_coord

Hi there,
I need to consolidate my local coordinates gotten in 10x10m sub-plots and to transform them into 100x100m reference. I am doing this by hand in excel... too much time!!
I need to sum up a reference measure to subplots measure by creating a loop..
The file "Coord_ha_R" has the columns names "Parcela Sub_parc X_local Y_local"
The reference file (ref_sub_parc_coord) to sum up values to X and Y_ local has  the colums name  "sub-parc X Y".

For each "XY_local fields in the "coord_ha_R", I need to sum up the respective XY  values in the  "ref_sub_parc_coord" file, considering each sub-parcel.

I'll really apreciate any help, because I have too many plots!!

Thanks
--
Kátia Emídio da Silva DSc
Eng. Florestal
Manaus/AM



Forestry Engineer
Manaus/AM-Brazil
Parcela	Sub_parc	X_local	Y_local
3	1	2.78	4.4
3	1	10	5.2
3	1	6.28	6.1
3	1	1.28	6
3	2	4.99	2
3	2	6.96	6.35
3	2	1.17	5.1
3	2	6.96	8.3
3	2	4.78	9.8
3	3	5.88	1
3	3	7.25	2.6
3	3	6.92	5.2
3	3	0.2	7.9
3	3	6.04	7.2
3	3	4.61	8.6
3	4	6.03	2.94
3	4	3.9	2.96
3	4	7.62	10
3	4	8.7	7.1
3	4	9.3	9.62
3	5	9.9	0.3
3	5	1.53	1.2
3	5	9.8	1.9
3	5	5.96	3.15
3	5	7.8	2.96
3	5	8.96	5.6
3	6	7.62	0.3
3	6	5.39	1.5
3	6	0.95	2.45
3	6	7.76	7.3
3	6	3.64	5.8
3	6	8.89	8.7
3	6	1.02	
3	6	5.31	8.5
3	6	2.09	8.6
3	6	0.46	9.3
3	7	0	3.55
3	7	3	5.8
3	7	10	8.05
3	7	0.73	8.8
3	8	2.63	0.85
3	8	4.2	4
3	8	5.37	4.5
3	8	0.96	6.7
3	8	7.03	8.1
3	8	7.32	9
3	8	8.6	9.7
3	9	6.91	1.3
3	9	4.53	3.3
3	9	7.32	7.3
3	9	1.09	6.1
3	9	0.54	6.4
3	10	5.54	3.3
3	10	10	8.2
3	10	5.32	5.7
sub-parc	X	Y
sub-parc1	0	0
sub-parc2	0	10
sub-parc3	0	20
sub-parc4	0	30
sub-parc5	0	40
sub-parc6	0	50
sub-parc7	0	60
sub-parc8	0	70
sub-parc9	0	80
sub-parc10	0	90
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo <at> r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Maurizio Marchi | 19 Jul 20:33 2014
Picon

PA selection in biomod2: how to merge different models before projection?

Dear List,
I have a question about biomod2 package.

I'm trying to follow instruction to select Pseudo-absences properly for
each Algorithm I want to use (GLM, GAM, MARS, MAXENT and RF) and I found
these publication very interesting:

Barbet-Massin M., Jiguet F., Albert C.H. & Thuiller W. (2012). Selecting
pseudo-absence for species distribution models: how, where and how
many? Methods
in Ecology and Evolution.

I was wondering: is there is way in biomod2 package to model each algorithm
separately (or some of them grouped together as suggested in the
publication) and then to merge them together to build an Ensemble model?
In the paper Authors suggested to average projections but I think that it
should be possible to do it with the BIOMOD_EnsembleModeling function

Many thanks, Maurizio

--

-- 
Maurizio Marchi, Ph.D. student
Florence, Italy
ID skype: maurizioxyz
Ubuntu 14.04 LTS
linux user 552742

	[[alternative HTML version deleted]]
Julie Lee-Yaw | 18 Jul 20:35 2014
Picon

"rasterize" shape file return empty raster

Hi

I have a shape file (initially produced by intersecting two other shapefiles using the PBSmapping
functions) and I would like to turn it into a raster based on a raster layer.

This has worked for me previously:

r<-rasterize(myShape, myRaster, mask=TRUE)

However, my end result is a "blank" raster. That is when I check values I get:

values      : NA, -Inf  (min, max)

Can anyone suggest issues that might cause the result to be a blank raster? I'm totally stumped given that
I was able to run this code for a different shapefile (also the intersection of two other shapefiles, one of
which is in common with those used to produce "myShape" above). I am also able to successful rasterize both
of the shapfiles used to produce "myShape" (myStates and myRange below) which makes me think it's a
problem with the intersection. To do the intersection, I used the following script:

library(PBSmapping)

#read in shapefiles
myStates<-reaShapeSpatial(file.choose())
myRange<-readShapeSpatial(file.choose())
p<- CRS("+proj=longlat +datum=WGS84")
proj4string(myStates)<-p
proj4string(myRange)<-p

#convert to PBS format
myStates.pbs<-combinePolys(SpatialPolygons2PolySet(myStates))
myRange.pbs<-combinePolys(SpatialPolygons2PolySet(myRange))

#intersection
myShape<-combinePolys(joinPolys(myRange.pbs,myStates.pbs,"INT"))

#convert to sp format for writing
myShape.sp<-PolySet2SpatialPolygons(myShape)

# and to spatial dataframe
IDs<-sapply(slot(myShape.sp,"polygons"),function(x) slot(s,"ID"))
df<-data.frame(rep(0,length(IDs)),row.names=IDs)
myShape.sp2<-SpatialPolygonsDataFrame(myShape.sp,df)

# write and then read back in before trying to rasterize
writeSpatialShape(myShape.sp2, "myShape")

myShape<-readShapePoly("myShape")

The intersection appears to work (e.g. plot properly) but one thing that I've noticed is that the joinPolys
command in PBS mapping no longer allows me to set maxVert...which may be more important for this
particular intersection...this is the only difference I can find and I'm not sure why it would matter.

Any help would be much much appreciated.

Thanks!

Julie

	[[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo <at> r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Alex Zvoleff | 18 Jul 20:18 2014

Unexpected behavior of raster "mask" function

I am using the mask function in the raster package (2.2-38) to mask out
areas within an image that are outside an area of interest (AOI). There are
NAs within this AOI that are meaningful - after masking the image I am
using freq to tabulate these NAs along with the other values in my AOI. For
this reason, I am using the "updatevalue" option in the mask function to
recode all areas outside of my AOI to a value (say 99) so that I can ignore
these areas in subsequent analysis, without having ignored areas share a
value (NA) with areas inside my AOI.

However, the mask function does not operate as I expected - areas in the
image that should be masked but that are NA prior to masking are not set to
"updatevalue". Is there an way to tell the mask function to recode all
areas (regardless of their initial value) to updatevalue? The below example
shows my problem:

library(raster)
# Setup an image with two NA values. Make one NA value inside
# the AOI, and make the other outside the AOI
img <- matrix(1, nrow=5, ncol=5)
img[3, 2] <- NA # Inside AOI
img[1, 2] <- NA # Outside AOI
img <- raster(img)
plot(img)

# Setup a mask where only the 3x3 region in the center of img is retained
msk <- matrix(1, nrow=5, ncol=5)
msk[1,] <- NA
msk[5,] <- NA
msk[,1] <- NA
msk[,5] <- NA
msk <- raster(msk)
plot(msk)

# Apply mask
masked_image <- mask(img, msk, updatevalue=2)

# Notice that cell (1, 2) of masked_image is NA, instead of the updatevalue
(2).
# This cell was not affected by the mask because it was NA prior to calling
mask.
# Is there a way to tell the mask function to set all masked areas to
updatevalue,
# regardless of whether a cell is NA prior to masking?
plot(masked_image)

Thanks,
Alex

--

-- 
Alex Zvoleff
Postdoctoral Associate
Tropical Ecology Assessment and Monitoring (TEAM) Network
Conservation International
2011 Crystal Dr. Suite 500, Arlington, Virginia 22202, USA
Tel: +1-703-341-2749, Fax: +1-703-979-0953, Skype: azvoleff
http://www.teamnetwork.org | http://www.conservation.org

	[[alternative HTML version deleted]]
Alessandra Carioli | 18 Jul 10:00 2014
Picon

semivariogram geoR

Dear bloggers,

I am using geoR package to do semivariograms (variog). I would like to specify 50km distances with a maximum
of 2000km distance and I am at loss on how to do it. The Coordinates of my shape file are projected longitude
and latitude

A         13.515094 	37.44574
B         8.664217 	44.82956
C         13.155760 	43.50457
D         7.393502 	45.73274
E         11.856243	43.53362
F         13.555405 	42.88887
G         8.187188 	44.87850 […]

Here is my code:

lat <- coordinates(shape.shp)[,1]
lon <- coordinates(shape.shp)[,2]
data <- cbind(lat,lon,fert)
dist <- dist(data[,1:2])
summary(dist)
breaks=seq(0,12,l=21)
v1 <- variog(coords = fert.l[,1:2], data = fert[,5], breaks = breaks,na.rm = TRUE)
v1.summary <- cbind(c(1:20), v1$v, v1$n)
colnames(v1.summary) <- c("lag", "semi-variance", "# of pairs")
plot(v1, type = "b", main = "Variogram") 

Any help on the matter is greatly appreciated.

Ale

	[[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo <at> r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Luis Miguel Royo Pérez | 18 Jul 09:19 2014
Picon

Re: R "frozen" in GRASS


On 16/07/14 18:56, Luis Miguel Royo Pérez wrote:
> Hi everyone,
>
> I'm new in this matter. I've intalled R GRASS and the packages needed 
> for use R and GRASS together. When I run GRASS I type this in the console:
>
> R --save
>
> And I get this:
>
> R version 3.1.1 (2014-07-10) -- "Sock it to Me"
> Copyright (C) 2014 The R Foundation for Statistical Computing
> Platform: i686-pc-linux-gnu (32-bit)
> R es un software libre y viene sin GARANTIA ALGUNA.
> Usted puede redistribuirlo bajo ciertas circunstancias.
> Escriba 'license()' o 'licence()' para detalles de distribucion.
> R es un proyecto colaborativo con muchos contribuyentes.
> Escriba 'demo()' para demostraciones, 'help()' para el sistema on-line 
> de ayuda,
> o 'help.start()' para abrir el sistema de ayuda HTML con su navegador.
> Escriba 'q()' para salir de R.
> Traceback (most recent call last):
>   File "/usr/lib/grass64/etc/wxpython/gui_core/goutput.py",
> line 759, in OnCmdOutput
>
> self.cmdOutput.AddTextWrapped(message, wrap = None)
>   File "/usr/lib/grass64/etc/wxpython/gui_core/goutput.py",
> line 1238, in AddTextWrapped
>
> txt = EncodeString(txt)
>   File "/usr/lib/grass64/etc/wxpython/core/gcmd.py", line
> 103, in EncodeString
>
> return string.encode(enc)
> UnicodeDecodeError
> :
> 'ascii' codec can't decode byte 0xe1 in position 39: ordinal
> not in range(128)
> Traceback (most recent call last):
>   File "/usr/lib/grass64/etc/wxpython/gui_core/goutput.py",
> line 759, in OnCmdOutput
>
> self.cmdOutput.AddTextWrapped(message, wrap = None)
>   File "/usr/lib/grass64/etc/wxpython/gui_core/goutput.py",
> line 1238, in AddTextWrapped
>
> txt = EncodeString(txt)
>   File "/usr/lib/grass64/etc/wxpython/core/gcmd.py", line
> 103, in EncodeString
>
> return string.encode(enc)
> UnicodeDecodeError
> :
> 'ascii' codec can't decode byte 0xf3 in position 25: ordinal
> not in range(128)
>
> Everything seems ok, but when I try to run these couple of commands 
> nothing happens:
>
> > library(spgrass6)
> > G <- gmeta6()
>
> I'm running GRASS 6.4.3 on a Ubuntu 12.04.
>
> Anyone could help me?
>
> Thanks per advance.
> -- 
>

--

-- 

	[[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo <at> r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Herry | 18 Jul 08:56 2014
Picon
Picon

image(rs) returns error after reclassify

Hiya,

I am getting the following error:

> image(rs)
Error in setValues(outras, m) : 
  length(values) is not equal to ncell(x), or to 1

rs was created from 
reclassify(rster, rcl, filename = 'landuse4BA', overwrite =TRUE,
progress='text')->rs
this finished without error. 

rster displays fine with image(rster)

any suggestions?
 - thanks 
Herry

--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/image-rs-returns-error-after-reclassify-tp7586756.html
Sent from the R-sig-geo mailing list archive at Nabble.com.

Gmane