Hot answers tagged

14

The area of a circular buffer is a monotonically-increasing function of buffer radius (on a planar coordinate system anyway). So a simple search strategy can find a radius R such that the area of the buffer of radius R clipped to polygonal region A is (up to some tolerance) s. The simplest search algorithm would just be a binary search. Start with two ...


6

You can do it like this: library(raster) # example data g <- getData('GADM', country='BRA', level=1) ext <- t(sapply(1:length(g), function(i) as.vector(extent(g[i,])))) colnames(ext) <- c('xmin', 'xmax', 'ymin', 'ymax') head(ext) # xmin xmax ymin ymax #[1,] -73.98971 -66.58875 -11.145161 -7.121320 #[2,] -38.23634 ...


4

One of many ways could be df$address <- with(df, paste(city, state)) df # city state address # 1 Lexington Kentucky Lexington Kentucky # 2 Cincinnati Ohio Cincinnati Ohio # 3 Indianapolis Indiana Indianapolis Indiana Or paste(df$city, df$state) instead of with(...). you need to have the full "address" for ...


3

Since the data is in a geographic projection, you cannot assign 30 to the resolution. If you look at the result of doing so, the r.ext raster has 1 row, 1 column and 1 value. So, that problem solved. Try defining the resolution in the call to raster. Since sometimes nodata begets nodata, I am assigning a value of 1 to all cells in the reference raster ...


3

Hm, I don't really see the problem. overlay works just fine on my machine. Here's the code. # (Q)=AnnPrec/[(Tmax+Tmin)*(Tmax-Tmin)]*1000 Q <- overlay(AnnPrec, Tmax, Tmin, fun = function(x, y, z) { x / ((y + z) * (y - z)) * 1000 }, filename = "q", overwrite = TRUE) And here's the resulting image including the range of values. range(Q[], na.rm = TRUE) ...


3

You can do this quite easily with base graphics, I don't see a real need for spplot or ggplot: Sample data: Spain and a random raster: require(raster) es = getData("GADM", country="ESP",level=1) es =es[-14,] # drop canary islands r = raster(extent(es),ncol=200,nrow=200) r[]=runif(200*200) Plot the raster and add the polygons - the alternative clips the ...


3

I would start here: library(raster) g <- getData('GADM', country='BRA', level=1) plot(g) You can extract the coordinates from g, but that is probably putting the horse behind the cart if you want to make a map. xy <- geom(g)


3

I want to add some details to Farid Cher's answer as this is a very common problem. Using amatch can do wonders, but with these Spatial objects you should not use base::merge and not access the @data slot. That would inevitably leads to a terrible mess (base::merge changes the order of records, and they would no longer match geometries). Instead, use the ...


3

I cannot recover your error but it seems that the parameters b=... l=... for band and layer are not present and I don't know what is in your data field a=class_raster, may be it is not a numeric type that fits into the raster band, that you want to fill with these data or NA could also be a problem. Here a little how to, I use to test the gdal_raster in ...


3

The problem is you are assigning metadata to the object, and probably even overwriting the existing metadata. This happens here: projection(pr) <- mycrs You should first run this to see what raster thinks it is already: projection(pr) or just print the object out to get a fully summary: pr Going out on a limb, I think you should warp this ...


2

Perhaps you need to first look at how ifelse works. I get the same results when I use it "stand-alone" and within a call to raster::overlay. a <- rep(2, 5) b <- rep(1, 5) d <- c(2, NA, 2, NA, 2) library(raster) r <- raster(nrow=1, ncol=5) A <- setValues(r, a) B <- setValues(r, b) D <- setValues(r, d) s <- stack(A,B,D) ifelse(a==2 ...


2

I downloaded your data and had a play around in QGIS - if I choose EPSG:102022 (Africa_Albers_EAC) it looks like it lines up. Depending on your plans for this that might be close enough - obviously for real science(tm) you'll need to go back to NOAA (or the publications listed) and check. Or you could try editing the spheroid definition below. ...


2

Ideally, when reducing by a factor, if there is a multimodal result I'd like aggregate to randomly assign the new cell one of the modal values and not always choose the same (if that's indeed what it does). That is not what it does. See ?modal and the ties argument. Your question is really about the modal function which you pass on to aggregate ...


2

if you need a separating character like a space or symbol (eg., Lexington - Kentucky ) you can use "paste" assigning the "sep" argument a value, otherwise "paste0" will join the strings with no seperator. df <- data.frame(city = c("Lexington", "Cincinnati", "Indianapolis"), state = c("Kentucky", "Ohio", "Indiana")) ( df <- ...


2

Define a grid at the desired resolution and rasterize by mean point value: library(raster) g <- raster(pts) ## gives an empty 10*10 grid on your extent/crs ## set the resolution (pixel size x/y) res(g) <- c(5, 7) ## whatever you want ## rasterize by "value" column with na.rm optionally r <- rasterize(pts, g, field = pts$value, fun = mean, na.rm ...


2

Include fill/color in your aesthetics mapping: ggplot() + geom_polygon(data = tz.prj, aes(long, lat, group = group), fill = "#f1f4c7", color = "#afb38a") + geom_point(data = tz.c, aes(lon, lat, fill = "Hospitals"), pch = 21, color = ...


2

I've just been doing the same thing. Pascal's answer almost covers it but you may need two extra steps as below. #After you create your list of latlongs you must set the proj4string to longlat proj4string(dat) <- CRS("+proj=longlat") #Before you re-set the proj4string to the one from sodo you must actually convert #your coordinates to the new ...


2

I am curious as to why you are not approaching this problem using a point pattern analysis? It is apparent that you are after a multiscale comparision but, it is not clear as to what end or what type of supported inference would be made. The type of standardization that your are attempting is hinting that a PPA would be a more supported methodology. ...


2

Please show your sessionInfo() and traceback() after the error occurs. This error message: Error in cbind(poly, rbind(poly[-1, ], poly[1, ])): trying to get slot "coords" from an object of a basic class ("NULL") with no slots is very odd as it would normally only occur when expecting SpatialPoints. As it seems you have very many, presumably very small, ...


2

Try clipping the polygons before using them (also, please try to provide complete code including library calls in the future): library(ggmap) library(rgdal) library(rgeos) library(ggplot2) URL <- ...


2

If I understand your specified model correctly, a separate covariate for each landcover type proportion. Yes, the spatial data structure will need to be addressed. You will need to create a separate raster for each landcover type that represents proportion of each given type at the appropriate scale. Depending on how your independent variables were derived, ...


2

I cannot recover your error, because do not know the content of the Corine_Hungary_2006.asc file. You could either work with the vector stuff you will find under http://www.eea.europa.eu/data-and-maps/data/clc-2006-vector-data-version-3 and select/ classify the right vectors together in a layer. Or you follow this how to: At first, get the corine data form ...


2

If you read your data into a raster package object you can use reclassify. Sample data, a grid of 1 to 25: > x=raster(ncol=5,nrow=5) > x[]=1:25 > as.matrix(x) [,1] [,2] [,3] [,4] [,5] [1,] 1 2 3 4 5 [2,] 6 7 8 9 10 [3,] 11 12 13 14 15 [4,] 16 17 18 19 20 [5,] 21 22 23 24 25 Then set ...


2

I'm not shure which verion you use but if you use a gdal version before 1.8, you have to create the tiff with a proper extention and resolution before you use it. http://www.gdal.org/gdal_rasterize.html ...look at dst_filename: The GDAL supported output file. Must support update mode access. Before GDAL 1.8.0, gdal_rasterize could not ...


1

This is how you can do that: library(raster) # erase x <- shp1 -shp2 # union (append) y <- x + shp2 But it won't be faster than gDifference (which is used under the hood --- the benefit of the method shown here is that attributes are not lost). Perhaps this speeds things up a bit agg <- aggregate(shp2) x <- shp1 -agg y <- x + shp2 You ...


1

This should work for removing NA's for a specific column yet, retain the sp class of the object. We will use the muese dataset from the sp library as an example library(sp) data(meuse) coordinates(meuse) <- ~x+y names(meuse) Here we add some NA's to copper, at rows 2, 5 and 20, and look at the resulting data meuse@data[c(2,5,20),]$copper <- NA ...


1

You just need to cut the data first: library(rgeos) library(maptools) library(raster) library(ggplot2) library(dplyr) library(ggthemes) library(ggalt) library(scales) #load shapefile ken <- getData("GADM", country = "KEN", level = 1) # make it less complex ken <- SpatialPolygonsDataFrame(gSimplify(ken, 0.001, TRUE), data=ken@data) #fortify for ...


1

It's almost impossible, due to the position of the points. You can create buffers of 400km2, but points closer to the coastline will always have a smaller area compared to the ones further away (>400km2). The only thing you can do is do perform a buffer analysis on the points and clip the created buffers with the coastline feature afterwards.


1

If you want to stick with spplot, you could e.g. use layer from latticeExtra to add 'Spatial*' objects to an existing plot. Based on the code provided by @Spacedman, this would look as follows. ## load package library(latticeExtra) ## plot raster and add polygons spplot(r, scales = list(draw = TRUE), col.regions = terrain.colors(100), at = seq(0, ...


1

You can pass a simple function, using "which" to return the desired julian day index, to the raster "calc" or "overlay" function. This will return a single raster layer with the first julian day of rain. He we create some example data that approximates your problem with 20 raster layers in the stack. library(raster) r <- raster( xmn=10, xmx=21, ymn=6, ...



Only top voted, non community-wiki answers of a minimum length are eligible