Dear Karl,
For the area of a spherical polygon, you could do
library(geosphere)
areaPolygon(xy)
I do not know how accurate this is relative to your alternative approaches.
Robert
Calculating areas on Earth’s surface
I’m interested in the distortion in apparent area from various
Pierre,
The trick is I'm using RColorBrewer, who does not provides this
posibility.
No, but the alpha argument to plot (with a Raster* object) does:
r - raster(system.file(external/test.grd, package=raster))
library(RColorBrewer)
col - brewer.pal(9, 'GnBu')
plot(r, col=col, alpha=0.3)
I'm trying to replace NA value with 0. But the function replace all
the value with zero:
dem[is.na(dem[])] - 0
Alfredo, that is indeed incorrect; fixed in version 1.9-5. Another way to do
this also works in 1.9-2:
dem[is.na(dem)] - 0
Thanks, Robert
--
View this message in context:
Apologies-- there was an error in the stack correlation function of my
last post. Below is an updated
version of the function.
Thanks Logan,
That is very useful, but I would not use stackApply, but the simpler calc.
Below a reworked version of your function:
stackcor - function(s1, s2,
, 2011 at 9:23 PM, Robert Hijmans [via R-sig-geo] [hidden
email] http://user/SendEmail.jtp?type=nodenode=6629146i=0 wrote:
Jan,
Any advice? Should I do this via Python or can this be done via R but
more efficient?
Interesting problem, that I had not noticed before. I think this is
because
Pierre,
There are two ways. You can use the alpha argument in plot, or in a function
that generates a vector of colors, see below. But note that this works in
the previous and in the next version of raster, but in the current version
only with 'oldstyle=TRUE'.
From ?slopeAspect
Jan,
Any advice? Should I do this via Python or can this be done via R but more
efficient?
Interesting problem, that I had not noticed before. I think this is because
a loop that can perhaps be omitted or else perhaps needs to move to a C
routine. I need to investigate that. For now, for
I need to crop a raster by pixel value, from upper X,Y left.
I'm using this:
spot - raster(.)
x_shift - 2
y_shift - 2
crop(spot,extent(xmin(spot)+res(spot)[1]*x_shift,xmax(spot),ymin(spot),ymax(spot)-res(spot)[2]*y_shift))
It works, but is there another way?
A little bit
I need to crop a raster by pixel value, from upper X,Y left.
I'm using this:
spot - raster(.)
x_shift - 2
y_shift - 2
crop(spot,extent(xmin(spot)+res(spot)[1]*x_shift,xmax(spot),ymin(spot),ymax(spot)-res(spot)[2]*y_shift))
It works, but is there another way?
A little bit
I have a process that creates a set of (large) raster layers.
Does it make sense defining an empty brick object
and then putting each layer as band of the brick?
Or is it better saving each layer and then make the brick?
Agus,
The simplest is certainly to write single layer raster files
When running the following script, I can plot the image.
Does anyone know how I can get the outline/border of the image?
Sample data is attached.
Muhammad,
I have no idea what you mean with how can I get the outline/border of the
image?, but perhaps the below helps:
library(raster)
r
I am trying to extract the per-pixel median from a raster stack, and
am finding that the process takes a very long time relative to taking
the per-pixel mean. It is long enough that I have never actually seen
whether it successfully completes on my data, but this dummy example
shows the
I'm trying to use subs function, but I've an error when I use by=
with the name of the column. But if I use the column number it work:
colnames(data_6SV_cluster)
[1] cluster intercept coefficient
subs(raster_cluster,data_6SV_cluster,by=cluster,which=2)
Error in `[.data.frame`(y, ,
I've a raster stack with 3 layers.
I need to classify it assigning one class to each unique combination of
values.
Perhaps the below works. It may be very slow if unq has many rows.
That was not quite right, I think it works when you do:
unq - unique(as.matrix(data_stack))
fun -
I simply would like to plot these on a map of the US. I've tried all sorts
of approaches, including
map('usa', project='albers', par=c(39,34))
points(x=data$LONGITUDE, y=data$LATITUTUDE, col='red')
You are plotting longitude/latitude data on an map with Albers projection.
Won't give you
Lines appear to be a GUI issue, went back to 2.13.0 and plot of raster
looks fine.
When plotting a Raster* object, raster calls the 'image' function. spplot
uses 'levelplot'. Both are affected by these white lines in R 2.13.1 on
Windows, I think (anyone seeing this on Mac or Linux? or not
I have a vector with repeated values like
...
I want to search the sequence in this vector. And search should give me
the
index of matching sequence.
If there is not a more direct way to do this, you can perhaps use the
movingFun function in raster for this.
library(raster)
values -
I have a jpg map of altitudinal contours in my study area which I
would like to map and use in R to calculate the altitude and slope of
GPS coordinates.
Estimating elevation from contours is not as common as it used to be because
of the availability of near-global high resolution raster
I am converting an .asc object into a raster object in order to use
functions from the raster package. I
then use the converted file to populate an array called FishModel for the
purpose of modeling. In the
process the columns and rows seem to be confused
Juliane,
It is easier to answer
I am using the dismo library for species distribution modelling using a
number of models. I would like to select different threshold levels. The
sdm.pdf vignette only gives an example to obtain the Maximum training
(or testing) sensitivity plus specificity threshold:
e@t[which.max(e@TPR +
I'm currently having some problems getting the crop() function in the
raster
package to work properly, and cannot for the life of me figure out what's
going on.
Essentially, I'm trying to crop a larger raster (raster1) to match one
with
a smaller extent (raster2), which is entirely
have a large number (several thousand) of asci files which were
delivered to me without a NODATA_value -999.0 line in the header
If all your files are for the same area you could do something like
r - raster(ncol=601, nrow=351, xmn=10, xmx=701000, ymax=25,
ymn=-101000)
filename -
Hi Eddie,
Consider this:
set.seed(0)
y = runif(365)
time = 1:365
m = lm(y~time)
m$coefficients
(Intercept) time
0.5297172419 -0.0001886475
But you do:
summary(m)$coefficients
Estimate Std. Error t value Pr(|t|)
(Intercept) 0.5297172419 0.0292052249
I am analysing monthly time series data using rasterbrick.
Commands below give me one slope plot
fun - function(x) { lm(x ~ time)$coefficients[2] }
slope - calc(s, fun)
plot(slope)
But what other 7 plots means when I run these lines...
fun - function(x) { lm(x ~ time)$coefficients }
I have a brick that contains 348 layers. Each layer is small in
dimension ncol=69, nrow=43. My question is, how do I calculate slopes
for each pixels from layer 1 to layer 348(Month one to months 348). I
am trying to find the trend as I am working with time series data.
for RasterBrick
From what I understood, most landscape classification methods available
from GRASS (through R) do not account for the temporaI aspect of
spectrum. So for example, if you have a raster where each band
represents NDVI values at a given date, from March to September for
example, the
Is there any way to
get a common color scheme
(i.e., setting common min and max values)?
Here is yet another way:
library(raster)
foo - raster(nrows=20, ncols=20)
values(foo) - runif(400)
bar - raster(nrows=20, ncols=20)
values(bar) - runif(400)
s - stack(foo, bar, foo*2)
# first
Interesting problem
You have to decide if you want to create 7 different shapefiles, one
for each of your grids and each containing one bounding box polygon,
or one shapefile with 7 features.
Here is a solution, I think, if you want a single shapefile with all
bounding boxes
I'm using the raster-package in R. I have problem writing an RasterStack to
a multi-band tif-file.
Error in attr(x@file, transient) - temp[[1]] :
no slot of name file for this object of class RasterStack
I believe this is a bug that has been fixed. Can you try again after
updating the
I have some rasters, respresentig a time series (let's say 5 years and 24
rasters for each year). The rasters are in a
raster stack or a SpatialGridDataframe, whatever is better.
Now I want to do some time series analysis which call for a ts class
object.
...
Now what I want to do is to do
I am working with crops planted area maps from two distinct sources.
One of the maps is based on a maximum NDVI composition, and the other
map uses joint information from satellite and census to estimate the
planted area.
Although the sources employ different methodologies to map the area
According to the help file, the predict() function in the raster package
creates a RasterLayer. Since a
RasterLayer cannot store more than one layer like a RasterStack, I wonder
if it is possible to use the
predict() functions when the output of the prediction is multivariate? For
example
Thanks for pointing this out. Some functions in 'raster' are indeed
slower/faster than others. Some seem to be reasonably quick. The focal*
functions are slow, and indeed, as Rainer alludes to, good candidates for a
C-implementation of at least the basic cases (e.g. fun=mean). Although I
recently
Thiago,
You can use the 'crop' function.
Best, Robert
Dear R colleagues,
Currently, I am working with the vegetation map published by SAGE
(http://www.sage.wisc.edu/download/potveg/global_potveg.html). This map is
published in netcdf and I can easily load it into R using the package
The second argument in projectRaster is another Raster* object, not a proj.4
string.
I expect this problem goes away if you do:
x - projectRaster(r, crs=+proj=longlat +ellps=WGS84 +datum=WGS84
+no_defs)
Robert
Having an error in getting library raster to read data from a netcdf
file when
I would like to overlay two raster objects to produce a nice map. The
underlying raster is a hillshade, the second is transparent.
plot(hillshade, col=grey(1:20/20))
plot(soilraster, alpha = 0.7, add = TRUE)
Does anybody know, how to hide the legend of the hillshade?
In addition to
I would like to confirm that the get.var.ncdf function within the ncdf
package automatically scales output data values if there is an embedded
scale factor within the ncdf file.
It does and that is easy to infer from the code frament below that you get
when typing get.var.ncdf. As you can
Lyndon,
I have no idea what might be causing this. I would suggest to track down
where exactly this happens, by going, line by line, through the code of
'subs'
getMethod(subs, c(RasterLayer, data.frame))
Presumably, something is not moving here:
for (i in 1:tr$n) {
v - getValues(x, row =
I am trying to resample a raster layer to another raster layer (dieffernt
resolution, different extent) using the resample function as implemented
in the raster package.
(...)
Interestingly, the same happens when using the example as provided in the
package description:
(...)
s - resample(r,
Does anyone know of a package (or a suggestion on how to implement) to
calculate, for two classified raster images of the same location but
different times, the relative probability of transitioning from one class
to
the other? Additionally, once this is figured out, how to apply this
Dear All,
I tried to write a function that can be used in focalFilter of raster
package, but the output
result is different from what was calculated using spdep for a sample
dataset.
It might be my mistake in the formula that is used in the function.
Following you will
find both
Dear Robert
I hope the following example is clear. you should be able to cut and
paste...
Nick
Dear Nick,
Thank you for the example. It reveals a serious error in focalFilter.
Results with that function in raster versions lower than 1.8-11 are wrong in
the first rows (half the number
Manuel,
bio1ai - crop(bio1, ai)
Error in intersectExtent(x, y) : Invalid extent
This somewhat cryptic error message means that the extent of y does not
overlap with the extent of x. Hence y cannot be used to crop x.
extent(bio1)
extent(ai)
Perhaps this happens because bio1 is in lon/lat
I am experiencing a strange problem when using the focalFilter
function in the raster package.
Nick,
Could you provide a self-contained example that does not require downloading
files? Given your description that should not be very hard.
Robert
--
View this message in context:
I have a raster file (cell value varies from 0 to 1) and a ESRI
shapefile (there are countries, and provinces inside each countries),
the two files overlapped in the same region. I want to summarized the
cells (with different value intervals, say 0.0-0.25, 0.25-0.5,
0.5-0.75, 0.75-0.1) by
Advait,
I encounter two errors; A) when using extent to crop the rasters
lims - extent(rcmraster,180,185,300,305)
Error in .local(x, ...) : unused argument(s) (180, 185, 300, 305)
I also tried using lat-long values in place of the index to no avail.
I cannot reproduce the problem.
I'd like to extract celle numbers from a raster using a line, but it seems
the cellnumbers=TRUE option isn't working with lines in raster::extract.
Etienne,
It is not working because cellnumbers is currently only an argument for
extracting values from a Raster* for polygons. That's what the
I tried your suggestion doing this:
library(raster)
myobject
-raster(C:\\Users\\komine\\Desktop\\MODIS\\MYD09GQ.A2010280.h16v07.005.2010282160347.sur_refl_b01_1.tif)
e - extent(-90, -32, -60, 15)
r -
Hi, I was wondering if someone could direct me to a package that would
allow you to query Google (or any other resource) for the latitude and
longitude of a ZIP code.
You can try the 'geocode' function in the dismo package.
geocode('03861')
ID lon lat lonminlonmax
I have an image which is a SpatialGridDatafFrame. This image has a number
of
cells that i need to erase or to delete but i don`t know how to do this in
R.
Mary,
It is not clear what you mean with erase or delete.
If you refer to setting some cell values to NA you can do, with
Rahul,
I wanted to extract the min value and max value from
the object. How can I
do this?
I have used the following function for this purpose
cellStats(r,'min')
[1] 0
But it is returning 0, not 0.001407638
The combination of
minValue(r)
[1] 0.001407638
I've got a SpatialGridDataFrame SG with density values given by spatstat
package
I want to get the average values of this grid for each polygon of a dataset
named iris
Mathieu,
I think you can do this:
library(raster)
#make a RasterLayer
r - raster(SG)
# extract mean values for polygons
Don,
The colortable is apparently not preserved when using crop() from the
raster package. This is something I would find useful.
I agree. Color table support is relatively recent and it is ignored by most
functions (and hence lost when using crop). That needs improvement.
(also, is there
I have a classified image of whole India (India measures 3,214 km (1,997
mi)
from north to south and 2,993 km (1,860 mi) from east to west) at 56m
resolution. I have generate a square grid of type polygon (area of each
square polygon in grid is 1sq km) that covers whole India. Now I have
I am getting the following error when I run aggregate (from the raster
package) on a raster image and output to a file.
Error in .GDALnodatavalue(dataformat) : cannot find matching nodata value
Ned,
I think this is something that has been fixed but had not made it to CRAN
yet. Could you
Im not sure what the problem is Robert, my previous post included both the
code file
and my albers coordinates. I will attach them again.
Got them now
Simulated data doesn't really prove anything to me if you can't apply the
results back to the earth.
? really?
This seems to work fine.
What are good raster formats and R libraries for dealing with non-square
cells?
You can use geotiff, as in my example (file.tif).
Is this proper use of the hpi values from R?
I do not know that
How can I effectively get these rasters back to square cells
I can see two ways to make the
Im having trouble with the results from the ks package being reversed.
There must be a simple error that I am overlooking? Does anybody see my
problem
Perhaps I could have seen what you are overlooking if you had provided a
simple self-contained example. That should be very easy to make in
Here is a self contained sample with a different set of points, same issues
though.
curset - read.delim(sample.csv, header = TRUE, sep = ,)
Thanks for trying, but this is not self-contained, as it refers to a file
that only you have.
When I use your routines they fail, saying that the
In order to illustrate the spatial domain of the General Circulation Model
(GCM) output data I am using, I
need to display some grids within a map. They need to be 2.5°x2.5° lat-lon
over specific coordinates.
Following, I create the empty grids to draw:
library(raster)
I discovered a strange behaviour of the function writeRaster in ther
raster package. Having a *.img or *.tiff file with a projected extend
object like this,
...
alway causes a wrong projection and associated coordinates when writing
to a envi file with the function call:
Carsten,
I'm trying to use MODIS 1day raster (.tif file) to obtain a, say, weekly
raster which is the average of the seven 1 day raster.
Using the following code (library raster) I can't obtain the results due to
computational limits.
Luca,
You can do something like:
s -
A small caveat - if the coordinates are planar, the explanation holds, if
geographical (degrees), it doesn't. The naive area values have been used
internally to plot larger Polygon objects before smaller Polygon objects,
to avoid overplotting in earlier versions of R that did not identify
I'm trying to read grib files using Raster, form the ECMWF interim
dataset (example file linked below). I don't get errors but the
coordinate system of the read raster does not seem to match.
Bart,
I think this is due to an instance of the grid vs. cell-center convention
mistake in gdal,
Steve,
But now I don' see an approach to subsetting this on the temporal domain.
Does raster support this kind of analysis?
x - subset(ss, 2:10)
x - ss[[2:10]]
But perhaps you do not need to subset, as you can probably use 'stackApply'
for your purposes
cellStats should be able to handle
I see functions to rasterize in the raster package but nothing to
vectorize...
and certainly nothing to do buffer around polygons.
raster has rasterToPolygons -- but it only works well for smallish rasters
(or when few cells have non NA values) because each cell with a value
becomes a
Brian,
I would like to know how to also reurn the ID of the polygon when using
the 'extract' function to extract raster cell values.
The values are returned in the order of the polygons. So the IDs would be
1:n, or in your example, poly$FID. Isn't that enough?
I am just extracting and
I was trying to do add the ability to ignore NAs (just like it's in the sum
method with na.rm = TRUE).
I can't find the way however.
You can also compare with the focal functions in raster; that should support
NA values
library(raster)
r = raster(meuse.grid, log.zinc)
rf = focal(r, 3)
Dear mailing list.. Please give me some information.
Which one is best fit for Remote Sensing Data RASTER or RGDAL package???
with Best Regards,
I would add that raster builds on rgdal, and so, if you use raster, in most
cases it would be raster rgdal, with rgdal providing the basic
Hi Peter,
I think you have made things way too complicated. I would suggest to do
something along these lines:
library(raster)
# first get the filenames.
vargroup - c(bio,prec,tmax,tmin,tmean)
filenames - list()
for (v in vargroup) {
filenames - c(filenames,
70 matches
Mail list logo
