Hot answers tagged

21

From the sp::over help: x = "SpatialPoints", y = "SpatialPolygons" returns a numeric vector of length equal to the number of points; the number is the index (number) of the polygon of ‘y’ in which a point falls; NA denotes the point does not fall in a polygon; if a point falls in multiple polygons, the last polygon is recorded....


20

You can get the same result by using st_join: First create a demo polygon and some points with sf. library(sf) library(magrittr) poly <- st_as_sfc(c("POLYGON((0 0 , 0 1 , 1 1 , 1 0, 0 0))")) %>% st_sf(ID = "poly1") pts <- st_as_sfc(c("POINT(0.5 0.5)", "POINT(0.6 0.6)", "POINT(3 3)")) %>% st_sf(ID =...


13

sp provides a shorter form to select features based on spatial intersection, following the OP example: pts[ply,] as of: points(pts[ply,], col = 'red') Behind the scenes this is short for pts[!is.na(over(pts, geometry(ply))),] The thing to note is that there is a geometry method that drops attributes: over changes behaviour if its second argument has ...


11

The simplest way to overlay two plots might be using the add = TRUE option in plot. Here is an example with artificial data # Load sp package for creating artificial data library(sp) # Create sample town points towns <- data.frame(lon = sample(100), lat = sample(100)) towns <- SpatialPoints(towns) # Create sample polygon grid grd <- GridTopology(...


10

Here is how I would do it in a Desktop. Get QGIS Desktop, an Open Source and popular GIS Desktop. Install Openlayers Plugin(In Menu, Plugins -> Fetch Python Plugins would display the plugin) Set Project CRS to EPSG:3857 Add Bing/Google Satellite Layer( Plugins -> OpenLayers -> Add Bing/Google) Add PostGis Layer(Layer -> Add Postgis Layer)


10

(Too long for a comment but thought this might be helpful in some way.) There is a plugin called Proportional circles which allows you create a legend broken into x number of segments (polygons) with a radius you can specify. You can download this from the menubar: Plugins > Manage and Install Plugins... Center your canvas to your points then click the ...


9

These are the basic steps you need to run through (is this what you have done? Apologies if this is too basic)... Make sure that your 1940 shapefile has an attribute containing county area (e.g. ‘Area1940’) use the Union tool (ArcToolbox > Analysis Tools > Overlay > Union) with your 1880 and 1940 shapefiles as inputs. Specify a new output shapefile. In the '...


8

PBSMapping should fit your needs. There's a tutorial at NCEAS. The code below is adapted from that tutorial. I'm making assumptions about your data btw. Please edit as appropriate for your situation. library(PBSmapping) #prepare towns pts <- read.csv("towns.csv") towns <- points(towns$lon, towns$lat) # read in shapefiles rivers <- importShapefile(...


8

I wanted to point out that you can rewrite the mean function, mean, which you can write yourself to do anything you want, including ditch the 0 values and calculate the mean. For example, if you want to ignore 0s: meanIgnoringZeroes <- function(x) { mean(x[x!=0],na.rm=T) } Then you can pass the function, meanIgnoringZeroes to overlay: mean <- ...


8

INPUT: After dissolving zones use following WORKFLOW: arcpy.Intersect_analysis("GRID #;ZONE #","D:/Scratch.gdb/intersect") arcpy.Sort_management("intersect", "D:/Scratch.gdb/sorted","Shape_Area DESCENDING") # DELETE MINORITIES USING GRID ID arcpy.DeleteIdentical_management("sorted", "ID") OUTPUT SHOWS "SORTED" AND GRID: Transfer dominant zone ID to grid ...


8

Think of the over function returning a query, now you have to do something with it. One easy approach is to pass an is.na() evaluation to a bracket index of the the SpatialPolygonsDataFrame object. This will return a subset of polygons that do not intersect the points. First, create some example data. require(sp) data(meuse) coordinates(meuse) = ~x+y sr1=...


7

It's not possible to mix geometry types in layers. You will need to add separate layers. Also, I'd recommend splitting the layers by topic, e.g. one layer with buildings, one with general zoning info. Don't just throw everything into one layer, otherwise you will have a mess to clean up later. Also styling is going to be easier with more layers. As @Willy ...


7

What you're trying to do is known as apportionment. This takes a numerical attribute value of a feature and divides it in some way between pieces that feature is split into. There are a number of different solutions. In fact, there's a specific tool for it in the Business Analyst extension, but you probably don't have that. Esri has also published a ...


6

Your line ras_sub0<-rasterize(subarea0,raster_bath) is just taking the index number of the polygons and assigning that to the values of the raster. If you want just the intersection of your polygon and the raster: subarea0_bathy <- intersect(raster_bath, subarea0) Update: As @GodinA notes, looks like intersect() sometimes doesn't return a raster ...


6

You were on the right track with over. The rownames of the returned object correspond to the row index of the points. You can implement your exact approach with just a few addition lines of code. library(sp) set.seed(357) pts <- data.frame(x=rnorm(100), y=rnorm(100), var1=runif(100), var2=sample(letters, 100, replace=TRUE)) ...


6

I could not reproduce your problem in R (similarly to jareks, in QGIS). My output map was ok, apparently. This is data from year 2000. This is your code (adapted) to reproduce all the three years (1990,1995,2000) of population density available on the very interesting link you provided. They seem all fine. library(raster) library(rgdal) italy_map <- ...


6

I think there is no better way than to handle zoomend event. To group markers by zoom levels, use L.LayerGroup: turning them on and off will be easier. Calling map.removeLayer() or map.addLayer() twice won't produce any errors or add a layer twice: there is an internal hash that prevents such things. So you can just have a bunch of if (zoom > 6 &&...


6

You cannot have R objects called "2000", so presumably these are fake names? Your example actually should work, so you may want to double check why you think that the results are incorrect. @aaryno's approach should work. I would do this: library(raster) s <- stack(r2000, r2001, r2002, r2003, r2004, r2005) x <- reclassify(s, cbind(0, NA)) r <- ...


6

Take a look at the raster function in the raster package. It will let you create a raster with a specified extent, number of rows/columns and resolution. Here I will use characteristics of your data summary to create a 100x100 raster within the specified extent. I am passing an extent object to define the x and y limits. You can also use the specific ...


6

Step 1 Make bit rasters for each of the unique classes. This can be a 1-band rasters for each class, or a single raster with a band for each class (e.g. GeoTIFF). If using GTiff, you can use the creation option NBITS=1 to conserve space. You may also want to consider twobit rasters to store three-valued logic where the third (e.g. 2) is NODATA, which would ...


6

Ok I think I follow the question. I attempted this with a test data set. I'm using a geodatabase so that the area is calculated (recommended). If you must use a shapefile, calculate a field with geometry for the shape_area before you do step 2. I have a polygon layer named Poly and a fishnet grid named FNET. Poly has a field (MTYPE) and values of 1 (maroon), ...


5

If your diagrams represent values from polygons, here's what I do: Create a new shapefile of centroids ("Vector->Geoprocessing Tools->Polygon Centroids"). This is a point layer with all the data columns from the original polygon. Now edit this point layer, and move the points to wherever you want each chart. And finally recreate the charts from this point ...


5

What you are describing can be accomplished using Extract by Mask. For example, the image on the left shows raster imagery with a polygon overlay. The image to the right shows the result of the extract by mask operation. By default, ArcGIS assigns NoData pixels as no-color. For display purposes, I reassigned the NoData pixel color from no-color to black ...


5

As an alternative to the first answer, this approach uses flag fields on the parcel layer to indicate if the parcel contains upland forest or wetlands: On the parcel layer, add integer fields "has_forest" and "has_wetland" Select (by attribute) the "upland forest" landuse features Select (by location) parcel features that intersect the selected landuse ...


5

I would recommend clipping the raster to the shapefile, then in the resulting raster you can look at the number of cells present for each of your classes. The area can be calculated by multiplying the number of cells by the area covered each pixel (cell size squared). It's a different approach than the equally valid solution offered above but from a ...


5

You'll want to to ensure that both of the layers are in the same Coordinate System. Add a new field in the BG Data, and calculate the area, bg_area. Run Intersect on your custom polygons and the BG data to create an intersection layer. Ensure that you have the area of your custom polygons calculated. Arc should automatically calculate this in Shape_Area. ...


5

The answer to this kind of issue should be to rasterize the information on to tiles on the server and simply serve the user tiles of the routes that they want to see. There are a plethora of options for doing this, but I think this train of thought is the best solution. Render the desired view on the server and only pass the tiles to the user, wash and ...


5

I got the tiles to overlay correctly. The problem was in the re-projection done by both ArcMap and QGIS. When I was checking the reprojected shapefiles in ArcMap and QGIS, they were overlaying correctly and had the correct SRIDs. So I imported the shapefiles in WGS84 in PostgreSQL using the SRID4326 with shp2pgsql then used ST_Transfrom to reproject the ...


5

Such a slight systematic shift is usually due to a lack of datum transformation before reprojecting the data. You should test the different transformation and your data will overlap correctly. I can't tell which one is best for you based on the information provided, but you can test it relatively fast. EDIT: if this doesn't work, you have two solutions: ...


5

Try the commented and reproducible example below! I use the sf package because is much faster than rgdal to open Shapefiles. Also, I used the great package mapview for interactive visualization in R. I transformed the points object to polygon's default projection before doing the intersection. Note: I downloaded the Points.csv file from Dropbox link first ...


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