Tag Info

Hot answers tagged

5

Use the rgdal package. If the shapefiles have a projection defined rgdal::readOGR() is recommended. This package provides bindings to the GDAL library (Geospatial Data Abstraction Library) and access to projection/transformation operations from the PROJ.4 library. This supersedes the shapefile read/write functionality in maptools.


4

This method uses the intersect() function from the raster package. The example data I've used aren't ideal (for one thing they're in unprojected coordinates), but I think it gets the idea across. library(sp) library(raster) library(rgdal) library(rgeos) library(maptools) # Example data from raster package p1 <- shapefile(system.file("external/lux.shp", ...


3

you need to actually reproject the raster into a geographic (decimal degrees) projection using "projectRaster" or "spTransform". Also look at CRS sp definitions that specify your desired projection string. The example in the help for the "projectRaster" is quite clear in how to do this. If you coerce your raster data into a SpatialPointsDataFrame object ...


3

Here's a simple reproducible example. I give 2 approaches, one using mosaic/merge and another that just does the initial rasterize at the total extent of all shapefiles combined. The results are the same. Added based on Jeffrey Evans comment: You also need to consider what to do when you have overlapping polygons. If you want to apply a function (e.g. sum, ...


3

This is not really a raster in longitude latitude, it's just arrays of values (including longitude and latitude). You can deal with these explicitly like this: f <- "F18-SSUSI_EDR-NIGHT-DISK_DD.20150107_SN.26920-00_DF.NC" library(raster) ## treat these not as rasters, but as arrays of values ## though raster() is extremely helpful in simplifying the ...


3

Here is a suggestion using ggplot. It can be improved, but it gives you the basic idea. library(maps) library(ggplot2) us.map <- map_data('state') # add PADD zones us.map$PADD[us.map$region %in% c("maine", "vermont", "new hampshire", "massachusetts", "connecticut", "rhode island", "new york", "pennsylvania", "new jersey", ...


3

Here you go. A couple of utility functions and then the meat in one function (and no for loops :)) islines <- function(g1, g2){ ## return TRUE if geometries intersect as lines, not points inherits(gIntersection(g1,g2),"SpatialLines") } sections <- function(sl){ ## union and merge and disaggregate to make a ## set of non-overlapping ...


3

Thank you for clarifying your question as it was previously quite unclear. You can read a multiband raster using the stack or brick function in the raster package and assign the associated RGB values to an sp SpatialPointsDataFrame object using extract, also from raster. Coercion of the data.frame object (which results from read.csv) to an sp point ...


3

You can use apply, which is actually the basis of the rowMeans function. If you are concerned that your row means are not correct because of NA's, just use the na.rm = TRUE argument in rowMeans. library (raster) r <- raster(nrows=20, ncols=10) r[] <- runif(ncell(r)) r[sample(1:ncell(r),10)] <- NA ( r <- as.matrix(r) ) # Count number ...


2

An alternative to render LiDAR data and RGB values in 3D is FugroViewer. Below, there is an example with sample data they provide. I used the file entitled Bmore_XYZIRGB.xyz which looks like this: When opening in Fugro Viewer select the corresponding fields available within the file (in this case, a .xyz file): Then, color the points using the RGB ...


2

A group aesthetic is missing: ggplot() + geom_polygon(data=fortify(regions), aes(long, lat, group=group)) Otherwise the last point of a polygon is connected with the first point of the next polygon. See also here: SpatialPolygonDataFrame plotting using ggplot Remove connecting lines in ggplot2 geom_polygon


2

Without looking too deeply into what you are doing, returning the number of iterations is fairly straight forward. I have now made it so that infinite loops don't happen and made the threshold one of the function parameters. hindex <- function(x, threshold){ iCounter <- 0 # Make counter 0 repeat{ head <- x[x>mean(x)] tail <- ...


2

We have one of these files on hand: ## dp is the root to our local data repository f <- file.path(dp, "data", "ftp.cdc.noaa.gov/Datasets/ncep.reanalysis2.derived/pressure/air.mon.mean.nc") library(raster) b <- brick(f, level = 1) b On my system that uses the "ncdf4" package, but it could also use "ncdf". The global Pacific-view extent looks ...


2

The colors are stored with the TIFF file, but either the legend class names are not in there or else they are not read by raster. Therefore you will have to manually copy the legend names from the Web site and manually create an informative legend, but you can automatically display the codes for the colors. The first color in the color table is the ...


2

You have to provide a color palette: library("raster") mapd <- raster("MCD12C1_T1_2011-01-01_rgb_1440x720.TIFF") plot(mapd, col=terrain.colors(255), axes=FALSE) The rasterVis package has lattice based methods for plotting Raster objects (spplot(...)): library("rasterVis") spplot(mapd)


2

It is fairly straight forward to set up the looping logic with an i,j index. However, I do not quite get your logic. What happens after (n - 4)? You can only calculate the adjusted mean to day 361. Why does calc or overlay, with movingFun, not work for you? That aside, addressing your question, the missing piece is that you index rasters in a stack using ...


1

Here is a code, a bit brute force, but your rasters do not seem to be that big. I give the code with some matrices as an example that you can run and see it works. Some modifications for your code are needed. Basically, the trick here is to use grep() to find the number in the raster names' vector and to subset the name of th variable as a string. Then ...


1

Here's a potential solution using R. Its runs quickly for a 1000x1000 raster. Not sure how it would scale to your 800MB file. library(raster) # Create a 1000x1000 pixel raster of land use classes 1-4 with some spatial structure rmeter <- raster(nrow=1000, ncol=1000, xmn=0, xmx=1000, ymn=0, ymx=1000) cosxy <- function(xy, T, A) { a <- 2 * pi / ...


1

That's not really resampling. Resampling would be converting all pixels in each grid square to whichever stat value (mean or max or whatever COUNTS not values since it's categorical data). But you want four values per cell, not one. It's categorical, so zonal stats won't work because it doesn't offer a count stat. I believe you'll want Zonal Histogram. Your ...


1

You can use the raster package: library(raster) dem = raster.open("a_dem') # plot with R library(rasterVis) plot3(dem) # convert the raster to a matrix (as Volcano) demmat = as.matrix(dem) Now, you can apply the same procedure as in r2stl/demo/Maungawhau.R z = demmat x = 1:dim(demmat)[1] y = 1:dim(demmat)[2] library(r2stl) r2stl(x, y, z, ...


1

Using my base raster: base<- raster("test.tif") plot(base) I got: For your breakpoints: breakpoints<-c(0,5,10,15,20,25) plot(base, col=terrain.colors(6), breaks=breakpoints) With help(plot), raster option, you can see examples of other descriptions for this R command. Another plot option: plot(base, col=colorRampPalette(c("red", "orange", ...


1

Unfortunately, the data has no projection, because it is satellite raw data captured during the flight around the earth. The first two subdatasets contain the lonlat values, but these are not understood by GDAL. So I tried to make a scatterplot out of the data: I extracted the first three datasets to XYZ format: gdal_translate -of XYZ ...


1

The projection you are setting (proj4string(pr) <- "+proj=longlat +datum=WGS84") is wrong. That should be the projection of the original data, not the target projection. You say that the data is NOT on a lat-lon grid, so you have to find out which CRS your data has. Looking at the .nc file, it's not contained there, so you have to look through the docs ...


1

After a little editing i used this script and its working fine. Now i just have to do it for each ID and for each day :-) library(SDMTools) library(rgeos) library(maptools) library(lubridate) #points <- readShapePoints("C:/Users/mlra/Desktop/FunWithDistance/DummyBird.shp") points <- tellus_2012_2013 lon <- as.vector(points$long) lat <- ...


1

It all depends on the level of accuracy that you need. A coarse approach would be the spherical law of cosines. This has issues with small distances - some say that it is around 1km, others say down to a few meters. A better approach would be the Haversine Formula. This works well, however, it doesn't take into account that the earth is not actually a ...


1

I think you are looking for the Intersection tool available in QGIS. In the following picture each feature of the blue polygon has a soil attribute and while the red polygon is the polygon of your clip. Using the Intersection tool in QGIS (Vector -> Geoprocessing -> Intersection) or in Processing (open the Processing toolbox and type Intersection) the ...


1

It seems to be a question not so related to GIS, that is since you are treating the raster as a "list" (or a vector) of values. Anyhow you can use the following workflow in R: Assume you have raster objects called r1 and r2. >r1Vals<-getValues(r1);r2Vals<-getValues(r2) # Extract values of r1 and r2 cells to verctors ...


1

As the others have mentioned, the "best" method is probably up to personal opinion. I'll make a stand for QGIS =) @Thomas has pretty much summed up what ArcGIS can do, QGIS is also capable of doing similar automated processing, exporting models to python scripts and offers a list of tools from a range of sources such as GRASS, SAGA and functions specific to ...


1

ArcGIS and Python go very well together. ArcGIS offers you arcpy, a python interface to all ArcGIS tools and more. Even if you are a beginner it is quite easier to learn. You can for example start building a model in Model builder and export the content as Python script. Afterwards you adoped the script the way you want (eg add loops, or other custom code) ...



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