Tag Info

New answers tagged

0

You will need to use the ESRI tool structure and include the files in your script folder if you would like to reference them directly by name like you are doing in your code. http://resources.esri.com/help/9.3/ArcGISDesktop/com/Gp_ToolRef/sharing_tools_and_toolboxes/a_structure_for_sharing_tools.htm Otherwise, you will need to append the location of the ...


0

A resolution of 1 is clearly not useful in this case. What you are seeing is a map of one raster cell (~ 9.5 - 10.5; 46.3 - 47.3) that has a single value (36). Had you provided more information, e.g. show(rMerge), it would have been very easy to spot. Do change the resolution and the results really should be different. But the resolution should be higher, ...


0

You want to use "merge" with the all = TRUE argument. However, be very careful if you want to keep the relationship between the data slot and the rest of the slots in the spatial object. You could take an approach where the external data is not actually joined to the spatial objects data slot but queried during whatever analysis is being applied. This is ...


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

I have made a rough code example, that uses simple structures and attempts to make the script easy to understand and follow. It is likely that it is inefficient and could be structured much better. The fact that you have 3 years is the key that we need to consider here (I will completely disregard the potential of a leap-year). Adding additional years is ...


0

The method suggested by rcs does what you want, but I thought I'd add to it to show how to join the counts back to the polygons. My interpretation is that you want to count the number of points within regions, where regions is defined by the NAME attribute of your shapefile. However, your shapefile has many features for each region (i.e. NAME), which is ...


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


1

The issue is that R is installed under Program Files by default, which is write-protected. Install R into another folder such as C:\R\R-3.1.3 and the problem goes away when you set it as the R path in QGIS. The solution is described here: http://wiki.awf.forst.uni-goettingen.de/wiki/index.php/R_installation This worked running QGIS 2.8.1 Wien on Windows 7 ...


0

The issue is that R installations are write-protected by default. Install R into another folder such as C:\R\R-3.1.3 and the problem goes away when you set it as the R path in QGIS. The solution is described here: http://wiki.awf.forst.uni-goettingen.de/wiki/index.php/R_installation This worked running QGIS 2.8.1 Wien on Windows 7 64-bit. I imagine it will ...


0

The issue is that R is installed under Program Files by default, which is write-protected. Install R into another folder such as C:\R\R-3.1.3 and the problem goes away when you set it as the R path in QGIS. The solution is described here: http://wiki.awf.forst.uni-goettingen.de/wiki/index.php/R_installation This worked running QGIS 2.8.1 Wien on Windows 7 ...


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

Get the coordinates of the cell centres and create a Spatial object: spts <- rasterToPoints(r, spatial = TRUE) Transform the points to your desired target: library(rgdal) llprj <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0" llpts <- spTransform(spts, CRS(llprj)) The values are already copieds as columns on this ...


1

If you are just after a metric of performance, this is a fairly straight forward type of analysis to specify. 1) specify a model, 2) predict the model(s) at the data, 3) apply an accuracy/performance metric based on observed vs. predicted. We can step through the process thus (note; this is a dummy example so, ignore the REML errors): First lets specify ...


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


0

It appears that you have Projected Coordinates there (not Latitude / Longitude aka GCS Coordinates). It probably wasn't clear to you that that was the problem. See this post. Converting geographic coordinate system in R


0

I' think it is not a matter of R ..only. Here a small sketch howto realize the referencing stuff under Linux. You can use GDAL to create a georeferenced axis parallel image from your heli position data, assuming, that your image has the dimension width and height and the bounding box is given by may be GPS-coodinates stored in UTM-33 shown in this table: ...


0

So I went for the library(sp). I now made a data frame with 0 or 1 if the line was contained in the polygon. I later use this data frame to subset the correct lines for each polygon for area calculation and comparison. for (r in 1:length(LALarea.x) ) { for (q in 1:length(LDdia.x) ) { ...


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


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


1

You need to change the R folder path to something like this: C:\Program Files\R\R-3.1.3\ instead of pointing to the bin folder. I just checked it on QGIS 2.8.1 and it works. Cheers.


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


0

In R, use hist3d from the rgl package. demo(hist3d) will show you a very low resolution example, commented on here: http://www2.warwick.ac.uk/fac/sci/moac/people/students/peter_cock/r/3d_bar_chart


0

What do the "coast", "ro" and "bc" objects have to do with your problem? The issue may lie in that you are using "readShapeSpatial". Have you tried readOGR in rgdal? If you are reading a polygon shapefile, readOGR will result in a SpatialPolygonsDataFrame object. If you in fact, do have a SpatialPolygons object and want to coerce into ...


0

It is quite simple: library("rgdal") polygons <- readOGR('path_to/file.shp', layer = 'file') class(polygons) >[1] "SpatialPolygonsDataFrame" >attr(,"package") >[1] "sp" poly_df <- as.data.frame(polygons) # do some staff with "poly_df" that doesn't support SpatialPolygonsDataFrame # then convert it to SPDF back again ...


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


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


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


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

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


0

Using the TIGER file from here: http://www2.census.gov/geo/tiger/TIGER2014/COUNTY/tl_2014_us_county.zip you can do this: # use rgdal to read read in the shapefile to preserve projection information library(rgdal) soco <- readOGR(".", "tl_2014_us_county") # use soco$STATEFP to find states, SpatialPolygonsDataFrame behaves like a data frame ne <- ...


-1

Check geosphere package distance function or fossil deg.dist function. You have data in degrees and need to translate it into meters or feet before doing clustering.


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


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


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


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


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


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


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


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.


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


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


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


0

This is exactly what the function raster::mask(x, mask) is for. It sets cells in x to NA when the corresponding cell in mask is NA. library(sp) library(raster) # Create some sample data r1 <- raster(nrows=40, ncols=40, xmn=0, xmx=2, ymn=0, ymx=2) r1[] <- seq(1, 100, length.out=ncell(r1)) r2 <- raster(outer(1:20,20:1), xmn=0, xmx=1, ymn=0, ymx=1) ...


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



Top 50 recent answers are included