New answers tagged

0

Don't know if I should have modified my previous answer. Perhaps, but that covers some things not in this answer, so I decided to leave it. Anyway, the code below works well for me. It looks for all of xmlNodes in the kml file that are called "Folder" and then sets the layer parameter of readOGR to that xmlValue. Tested on working directory with around 6 ...


0

Don't know if this is still a problem for anybody else, but I was running in circles for a while with this. What finally worked for me is below. It uses the XML package to get at the xmlValue of the right node. I had to set the layer parameter of readOGR to the name of the one of the folders within the kml file. When I set the layer parameter to the of ...


0

You could approach this as a raster/vector integration problem, where your course resolution data are essentially sampling blocks represented as polygons. First, let's create your example data. library(raster) library(sp) ext <- extent(0,1000,0,1000) r1 <- raster(nrows=1000, ncols=1000,ext) r1[] <- sample(seq(from = 1, to = 6, by = 1), size ...


1

I believe that you can specify an IDW with a radial basis function in the geosptdb library using the "rbfST" funciton. The package is intended for spatial-temporal IDW but should run with a single temporal dimension. This is at least a good place to start and you may be able to modify the "rbfST" function to suit your specific needs.


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, ...


0

Try this, use "raster::intersect" and use the "group" for the map_id for ggplot: library(rgdal) africa_countries <- readOGR(dsn="shapefiles", layer="AfricanCountries", stringsAsFactors=FALSE, verbose=FALSE) precip_contours <- readOGR(dsn="shapefiles", layer="precip_contours", ...


0

I was able to use the randomForest package to answer my question which I explain in more detail here: http://stackoverflow.com/questions/34864182/predict-estimate-values-using-randomforest-in-r


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, ...


1

For example: library(ggplot2) ggplot(dfexp, aes(lon, lat)) + borders() + xlim(c(144.325, 144.40)) + ylim(c(-41.575, -41.5)) + stat_density_2d(aes(fill = ..level..), geom="polygon") + geom_point(position="jitter", alpha=.2, colour="white") Or using ggmap (as requested): library(ggplot2) library(ggmap) map <- ...


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 ...


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 ...


0

For US-level analysis I usually use the USGS Albers Equal Area projection "USA_Contiguous_Albers_Equal_Area_Conic_USGS_version". proj4string: +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs EPSG: SR-ORG:7480


0

Here are two approaches to compute correlation coefficients with Raster objects (notwithstanding the comments on your question about utility; that I concur with). # generate some data library(raster) set.seed(89) b <- brick(system.file("external/rlogo.grd", package="raster")) d[] <- runif(ncell(b)*nlayers(b)) s <- stack(b[[1:2]], d[[1:2]], b[[3]], ...


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 ...


0

Here would be a procedure suggestion using Matlab and the mapping toolbox. this is not a turnkey solution, but if you are familiar with matlab, it would provide you the advantage of automation and replication on many rivers. 1- Produce a centerline, either by hand in ArcGis or produce it automatically (google for a method). I have'nt found a way that ...


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 <- ...


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) ...


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 ...


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. ...


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 ...


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, ...


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 ...


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, ...


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 = ...


0

One solution would be to convert each of your polygons to rasters. Then take the portion of cell counts in each landcover class by the total cell count from the clip. Here is an example function: PortionClassInPoly <- function(MySingleShape, MainRaster) { # Start with one shape from your list # Rasterize this shape according to the ...


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 ...


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 <- ...


0

I also made a quick easy and holistic tutorial you folks can check out here https://github.com/mattjbayly/MapsProj. All you need is the R text script and all other material is downloaded remotely from R. My tutorial covered GIS data imports/exports, basic manipulation of vector & raster data, some of basic plotting and a brief overview of projections. ...


0

Not extremely efficient, but gets the job done city <- c("Lexington", "Cincinnati", "Indianapolis") state <- c("Kentucky", "Ohio", "Indiana") df <- data.frame(city, state) df[,3] <- cbind(paste(df[,"city" ], df[,"state"], sep=" ")) colnames(df[,3]) <- "address" > df city state address 1 Lexington Kentucky ...


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 ...


1

Looked at your R solution... Guess I'm overly visual..like to see what I'm doing. If you mess around with Geostatistical Wizard in ArcGIS you can interactively see the effect of various interpolation choices. The density thing is always an issue for the bathymetry sounding data I get. My solution is to create a tin from the points and a shoreline, as well as ...


0

I tend to use gArea from rgeos for such kind of operations since I'm typically interested in a single value rather than an entire raster layer with area information. To get area information in meters, I would recommend projecting your raster (which is obviously located somewhere in France) to EPSG:32631 (UTM zone 31N) before transforming it to ...


1

Is R actually crashing or is the code just failing? These are very different outcomes. I believe that your problem was that your data is in a geographic and not planar projection. I reprojected the data, from your link, to the USGS definition of Albers and gBuffer ran with no problems. You will need rgdal along with sp for the projection transformation. ...


1

Use the ?area function. library(raster) r <- raster() ## boilerplate longlat raster ar <- area(r) ## applies a correction for latitude, to km2


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 ...


1

In QGIS, you can use expressions to add fields using the Field Calculator Following expressions will give you the min and max coordinates for the polygon xmin($geometry) xmax($geometry) ymin($geometry) ymax($geometry) For center coordinates, you can use xmin(centroid($geometry)) ymin(centroid($geometry))


0

I solved this problem with these instructions: myData$longitude = as.numeric(as.character(myData$longitude)) myData$Magnitude = as.numeric(as.character(myData$Magnitude)) myData$latitude = as.numeric(as.character(myData$latitude))


1

In a bimodal distribution the fun=modal argument will result in assigning the value with the peak frequency in the distribution. Not sure what other behavior you expect? Statistically speaking, I do not see any support for randomly selecting a different peak in the distribution just because the distribution is bimodal. The only way I can see justifying ...


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. ...


1

When asking a question please show some of the code you have been using. In fact, you should try to create a self-contained reproducible example (there are many examples on how to that on this website). Otherwise it is almost impossible to help you. I am guessing that you are using the rasterize function from the raster package. The grid cells are rather ...


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 ...


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 ...


0

If a loss in the precision of the overlay is not terribly important - assuming it is precise to begin with - one can typically achieve much greater zonal calculation speeds by first converting the polygons to a raster. The raster package contains the zonal() function, which should work well for the intended task. However, you can always subset two matrices ...


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 ...


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 ...


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)


0

Not sure if this qualifies as a "trick", but I'm a huge fan of the combination of the acs package (for selecting US Census data) and the leaflet package (for making interactive javascript maps that can be hosted online). This tutorial does an excellent job illustrating the benefit of using these two packages together.


0

I ended up hosting the output of htmlwidget::saveWidget on a Github Pages project website. I embedded the html files in the ioslides presentation using an iframe, as demonstrated below: <iframe title="My Map" width="980" height="400" src="http://username.github.io/reponame/maps/my_map.html" frameborder="0" allowfullscreen></iframe> This ...


5

In QGIS, you could use the Polygon from Layer Extent... tool from the toolbar (Vector > Research Tools > Polygon from Layer Extent...). This essentially outputs a bounding box layer for each feature (if you select the option) with fields containing the coordinates of the max, centre and min of X and Y along with few other statistics:



Top 50 recent answers are included