New answers tagged

1

This is because a hole should always belong to a non-hole - i.e. you cannot have a Polygons where the only ring is a hole. Note that the hole argument is not ignored by Polygon(), it is overridden and the coordinate order reversed by ?Polygons() - "In Polygons, if all of the member Polygon objects are holes, the largest by area will be converted to island ...


2

Here's a way, use the raster package. (If you really have data for this it's probably not a great idea to grid in longlat and then resample to Mollweide - but depends on your purpose). R library(raster) library(maptools) r <- raster(surface) data(wrld_simpl) prj <- "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84" w <- spTransform(wrld_simpl, CRS = ...


0

You can do the analysis described in your post without converting the raster to a polygon. Use the raster::extract function to extract the raster values to each polygon. You can then use lapply on the resulting list object with table to return cell counts of each class. For area of each raster class, you just use a standard conversion of cell area and ...


0

from : https://sites.google.com/site/eospansite/installinggdalonmacosx I'm copying these things in case link be broken: Installing GDAL on Mac OSX and Polygonize Rasters Go to http://www.kyngchaos.com/software/frameworks and install the Unix Compatibility Frameworks in this order: GDAL Complete GSL framework FreeType cairo GDAL Complete: you ...


3

The package PostGIStools can help. See for example the vignette. Another way could be to transform your Spatial*DataFrame geometry to WKT, insert into PostGIS using the classic RPostgreSQL package and re-create the geom there.


0

You can use GRASS in QGIS as described here. Just create a slope layer then, using the raster calculator, make a new layer of a range of acceptable slope angles. You can then run the r.cost and. r.drain tools as mentioned in the other post and then polygonise the r.drain output and you've got your road. Remember that when you run the r.drain tool that you ...


1

The answer to my own question is that the local firewall settings prevented me from downloading the tiles from OSM. The way how to solve it was to set proxy (from within R > link). The con of this solution is that the viewer still doesn`t deliver tiles as the proxy settings at the moment do not apply to RStudio Viewer. Although, the leaflet output appears ...


0

What is the hypothesis of spatial process that you are testing? I am dubious over the logic of your analysis. What you describe is inherently a meta-model problem where error and spurious spatial is propagated, in unknown forms, through each step. This makes it impossible to quantify any associated significance. If this raster is the result of a local ...


1

I can not help you within R, but using GDAL takes you further: gdalinfo acpcp.2000.nc tells you that the first two bands contain the lat and lon coordinates, and they are in WGS84 degrees: SUBDATASET_1_NAME=NETCDF:"acpcp.2000.nc":lat SUBDATASET_1_DESC=[277x349] latitude (32-bit floating-point) SUBDATASET_2_NAME=NETCDF:"acpcp.2000.nc":lon ...


0

due to the comments i did this. >cord<-projection(stk, asText = F) >cord CRS arguments: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 But now whe i used cordas an input in the function raster(nrows,ncols,CRS=cord) i got an error because cord had to be text, so i did this. >cord<-projection(stk, asText = T) and finally ...


1

The rioja package provides functionality for constrained hierarchical clustering. For what your are thinking of as "spatially constrained" your would specify your cuts based on distance whereas for "regionalization" you could use k nearest neighbors. I would highly recommend projecting your data so it is in a distance based coordinate system. require(sp) ...


0

Have a look at mclust package for model based clustering. It supports higher dimensions (more than 2). Here is an example: library(mclust) #unclustered plot(iris[,1:2]) mod1 = Mclust(iris[,1:2]) plot(mod1, what = "classification") 3 dimensional would be: mod1 = Mclust(iris[,1:3]) plot(mod1, what = "classification")


0

CRS is an Interface class. It is not supposed to work like that. Use projection(raster) or proj4string(raster) as per comments.


2

A simple function can do that for you. Here's how: allKmlLayers <- function(kmlfile){ lyr <- ogrListLayers(kmlfile) mykml <- list() for (i in 1:length(lyr)) { mykml[i] <- readOGR(kmlfile,lyr[i]) } names(mykml) <- lyr return(mykml) } use it with: kmlfile <- "se\\file.KML" mykml <- allKmlLayers(kmlfile)


0

Hourrah !! There is a function to create only OUTSIDE buffer !! my first (working) function in R !! :) # create basic function outerBuffer<-function(x){ buff<-buffer(x, width = 1, dissolve = T) e<-erase(buff,x) return(e) } # create data p1 = readWKT("POLYGON((2 2,-2 2,-2 -2,2 -2,2 2))") # apply function, create only outside buffer ...


5

There is no way to define "outside only" in gBuffer. You have to go through the additional step of turning the inner polygon to null, and for a good reason. You can use the raster::erase function to remove the internal polygon. If you really want this as part of the gBuffer function why not just write your own modification of gBuffer that adds an "outside ...


1

If the original data was rescaled to 8-bit it should be 0-255 and not 0-200. That aside you can take a normalization approach but shift the centrality over so the distribution will bound into the negative. Two normalization formulas that will do this are: Formula 1: [(x - "x min") / ("x max" - "x min") - 0.5) * 2] Formula 2: ["new min" + (x - "x min") * ...


0

It is very unclear as to what your problem is! In the future please make an attempt to provide a clear problem, question and what you have tried (along with code). Here is a very general introduction to list objects in R that will likely, indirectly, answer your question. So, you apparently have a vector as an element in a list object. You can examine the ...


0

You can rotate the coordinates of an SP object using elide from maptools, but you'll lose any true geographic location reference. You also won't be able to overlay on map tiles since they are fixed at that alignment. Example using scot_BNG from readOGR: > plot(scot_BNG) > plot(elide(scot_BNG, rotate=-45)) Alternatively use an oblique mercator ...


1

Your raster and your polygon do not overlap. If they do: > r1 = raster(extent(c(107,110,38,40))) > polygon_mus <- extent(c(107, 111, 37, 40)) > cell <- extract(r1, polygon_mus) Then that works. If they don't: > r1 = raster(extent(c(112,114,38,40))) > cell <- extract(r1, polygon_mus) Error in (function (classes, fdef, mtable) : ...


0

I used rasterToPolygons from the raster package too in the past, but now I prefer gdal_polygonizeR by John Baumgartner. It bases on gdal_polygonize.py and is much faster. John Baumgartner published the code and gave an example for usage in his blog. If you are familiar with python you could use gdal_polygonize.py directly of course.


2

Gdal Installation Install Gdal command line tools and check to see if its binaries are added to path environment variable. e.g. in windows: open Run and type: rundll32.exe sysdm.cpl,EditEnvironmentVariables Then follow the screenshot Download and install gdal python bindings from here according to your python and OS. install it using: pip.exe ...


0

There is a workflow to fill in NA values by the nearest non-NA value in R here. You can easily set a maximum distance (in map units/meters).


2

This may be easier using the Orfeo toolbox (https://www.orfeo-toolbox.org/), this is provided with OSgeo4W and can be accessed usign QGIS or a command line interface. This tutorial uses mean shift segmentaion to generate objects, which can be the classified using SVM/random forests etc. ...


5

You need to make a spatial points object and then reproject that: library(sp) library(maptools) data(wrld_simpl) wrld_simpl.trans <- spTransform(wrld_simpl, CRS = ("+proj=moll +lon_0=0 +x_0=0 +y_0=0")) plot(wrld_simpl.trans, col = 'lightgray')) latitudes <- c(-30, 0, 20, 50) longitudes <- c(-20, 10, 50, 80) Make an SPDF: pts <- ...


0

You have to treat the @data slot (the data frame) a bit differently to reduce the number of attributes / columns: Set the working directory and load the right library require(rgdal) setwd("~/workspace/TEMP/") Make a fake SPDF to work with, as OP describes pts <- data.frame(x=1:5, y=1:5, A = 1:5, B = 1, C = 2, D = 3) coordinates(pts) <- ~x+y ...


0

To find the batch_gdal_translate function, a solution is to use the SOS package directly in R/RSudio to find it. library(sos) unique( findFn("batch_gdal_translate")$Package ) found 2 matches Downloaded 2 links in 1 packages. [1] "gdalUtils" If you use directly findFn("batch_gdal_translate") You go directly to batch_gdal_translate {gdalUtils} from the ...


0

It seems the size data will be used to select the surrounding pixel and around the center pixel of your choice. Size = c(1,1) will give you the output not only of the center pixel(or the point) but also 8 neighboring pixels. So if you want to download the data only for the concerned point then give 'Size' value as c(0,0). See the Vignettes for more detail. ...


-1

I implement a method to split polygon along the first two points using R. See function way_points here: https://github.com/byzheng/missionplanner/blob/master/global.R


0

If you're willing to use TileMill and MBTiles, you can do this pretty easily using CartoCSS to style the original Tiff, and output an MBTiles file. #dem { raster-opacity:1; raster-scaling:lanczos; raster-colorizer-default-mode: linear; raster-colorizer-default-color: transparent; raster-colorizer-stops: stop(0,#aaf) stop(100, #afa); } ...


0

I think you need to subset your dataset using bracket notation. As an example: point_pattern <- point_pattern[point_pattern$A == 1,]


3

I emphasize that this function might not be useful for everyone. It will depend on what you're trying to do, and your inputs, and how possible it is to do what you're trying to do given your inputs. For instance, if you have a complicated shape with concave sides, forget about it. In my case, I just wanted to cut a box filled with random points down to a ...


0

Here is a reproducible example using sp internal dataset meuse since I don't have access to your shapefiles: r1 = cbind(c(180114, 180553, 181127, 181477, 181294, 181007, 180409, 180162, 180114), c(332349, 332057, 332342, 333250, 333558, 333676, 332618, 332413, 332349)) r2 = cbind(c(180042, 180545, 180553, ...


0

A possible reason might be that vi (which you pass on to whittaker.raster) does not exist. Please provide the error message in order to clarify this issue. Note that MODIS has been moved to GitHub a while ago. You can easily install the latest versions via (stable version) devtools::install_github("MatMatt/MODIS") (development version) ...


1

You just need to loop over the combinations you need. > combn(1:length(Sl),2) [,1] [,2] [,3] [1,] 1 1 2 [2,] 2 3 3 that gets you (in columns) the line indexes. So then do: MyLines = apply(combn(1:length(Sl),2),2, function(x){ gIntersection(Sl[x[1]], Sl[x[2]]) }) Then MyLines[[1]] is the ...


1

The output of extract is a list of cell numbers and values stored in matrices, it needs tidying to a data frame. (Data frames are best since they can store different types of data, like integer cell numbers and numeric values - which is essentially what cellnumbers=TRUE is for). library(raster) r <- raster(volcano) ## simplify the values r <- (r %/% ...


0

Solution: I changed the order of the coordinates in the netcdf file, using the NCO operators (http://nco.sourceforge.net/nco.html#ncpdq-netCDF-Permute-Dimensions-Quickly) ncpdq -a lat,lon in.nc out.nc I am guessing that the R raster package expects the latitude variable to come before the longitude variable, although I haven't run across that in the ...


0

changing the function in the script to MAX, has solved the problem : m1 <- mosaic(r1, r2, fun=max)


0

Ok I use this NDVI.slope[NDVI.slope>0.03] <- NA In stead of NDVI.slope_2<-trim((NDVI.slope), padding=0, values = c(0.02:0.03)) ................ and It work !!!!! ...................


0

Unfortunately, your test example is not online anymore, so I took a gobal Land-Sea-Temperature sample from ftp://landsaf.ipma.pt/pub/Example_of_Product/ As mentioned in my answer to Transforming geostationary satellite image to lon/lat, the size of global geos reprojected files is 3712x3712 pixels (3-km-resolution), but the actual data array is smaller. ...


0

I resolved it using @mdsumner suggestions plus some tweaks: As suggested by mdsumner I used calc to get which values where of the values I needed, however this only gave me a 0 or 1 value, and I needed a 0 to 1 continuous landcover area type of thing, so then I used aggregate, and just selected how coarser my raster would be. It worked just as I needed with ...


2

Try fs <- list.files(path="F:\\MODIS\\Modis EVI\\HDF8 EVI", pattern = "tif$", full.names = TRUE) library(raster) s <- raster::stack(fs) writeRaster(s, "hdf8_EVI.TIF") I have no idea what stackSave is. Please read about asking questions here and look into the basics of the raster package. Don't use assign, and don't use setwd - both are really bad ...


1

As of one week ago, this can now be achieved with the development version of mapview: # devtools::install_github("environmentalinformatics-marburg/mapview", ref = "develop") library(mapview) library(sp) # video pt <- data.frame(x = 174.764474, y = -36.877245) coordinates(pt) <- ~ x + y proj4string(pt) <- "+init=epsg:4326" mapview(pt, popup = ...


0

I recently had to do a similar map. The solution I came up with uses the rasterVis package, rather than ggplot2 (which is an awesome package, by the way). In my case, I had a map of trends over time (which is a numerical variable too) and also a map showing the significance of the trend (obtained from a statistical test). In my case, I wanted to plot the ...


3

"I couldn't find it in any of them" ? CRAN Task View: Analysis of Spatial Data (Geostatistics) You can also use the SOS package directly in R library(sos) # find packages with a "spatial regularization" function unique( findFn("spatial regularization")$Package ) found 43 matches; retrieving 3 pages 2 3 Downloaded 29 links in 11 packages. [1] ...


0

I think the goodness of fit measure you are looking for is a chi-squared. For which (I believe) your 'minimum weighted sum of squares' IS your goodness of fit statistic. Assuming this or this is what you are doing. Reading about the packages on CRAN can sometimes be helpful. This thread will also be relevant. Or you could google: "calculate R2 weighted sum ...


3

Ok, this is not really an answer to the my question, but if you are interested in the biomass map I found a link from where you can download the actual data: http://whrc.org/publications-data/datasets/pantropical-national-level-carbon-stock/


1

Here's a potential approach: I haven't tried this so your mileage may vary, but it should be quite doable (but it might require rather heavy-duty scripting). I'm starting to think in a similar direction for some projects I'm working on. Ingredients: QGIS PostGIS (likely the best storage medium for relational data but others could potentially be used) ...


3

UPDATE: The above-mentioned approach by Oscar Perpiñán found here seems to work, if I only take a subset of the matrix the way Oscar did. For that to work I needed to know the subset boundaries and the static files, both available online for a few regions (such as North Africa). m = m[700:1850, 1240:3450, drop = FALSE] lon = raster("NAFR_LON.img") lat ...


7

Your attempt is designed to fail. If you look at the image, you see the data arranged as a circle, with black triangles in the corners of the square, where the satellite view goes right into orbit. In your test data, you see only NODATA -32768 for those parts of the image. The extent is between +/-75 and +/- 78, but these values are only reached in the ...



Top 50 recent answers are included