The Stack Overflow podcast is back! Listen to an interview with our new CEO.

Hot answers tagged

38

I would drop using the maps package and find a state shapefile. Then load that into R using rgdal, and then do some polygon overlay work. library(raster) # use state bounds from gadm website: # us = shapefile("USA_adm1.shp") us <- getData("GADM", country="USA", level=1) # extract states (need to uppercase everything) nestates <- c("Maine", "Vermont", ...


38

1) read your shapefile with Fiona, PyShp, ogr or ...using the geo_interface protocol (GeoJSON): with Fiona import fiona shape = fiona.open("my_shapefile.shp") print shape.schema {'geometry': 'LineString', 'properties': OrderedDict([(u'FID', 'float:11')])} #first feature of the shapefile first = shape.next() print first # (GeoJSON format) {'geometry': {'...


32

Here is an approach using extract() from the raster package. I tested it with altitude and mean temperature data from the WorldClim website (I limit this example to altitude, temperature works similar), and an appropriate shapefile of the US containing state borders is to be found here. Just download the .zip data and decompress it to your working directory. ...


29

The Raster|Extraction|Clipper tool will help you to do this. You can open the tool and then click and drag in the raster image to select the area you want to export as a new raster (Clipping mode: Extent), and then refine the exact coordinates in the Extent fields (if necessary). Probably no way to do contours on only a selection of a raster; clip it this ...


27

Use the Erase (Analysis) Tool:


23

I have finally gotten around to improving this function. I found that for my purposes, it was fastest to rasterize() the polygon first and use getValues() instead of extract(). The rasterizing isn't much faster than the original code for tabulating raster values in small polygons, but it shines when it came to large polygon areas that had large rasters to be ...


22

There is a QGIS plugin called Digitizing Tools: The documentation says: Split selected features with selected line(s) from another layer applies to: line and polygon layer (multi or single part) Splits all selected features of the active layer with the selected line features of another layer. The splitting creates new features (not multi features). Each ...


22

Since today, there is a st_crop function in the github version of sf (devtools::install_github("r-spatial/sf"), probably on CRAN in the near future too). Just issue: st_crop(nc, c(xmin=-82, xmax=-80, ymin=35, ymax=36)) The vector must be named with xmin xmax ymin ymax (in whichever order). You can also use any object that can be read by st_bbox as ...


21

You need to loop through your inputs. Multivalue is semicolon delimited. Split on that and loop through them. (AddMessages to show how the fcs are presented) import arcpy ins = arcpy.GetParameterAsText(0) arcpy.AddMessage(ins) for fc in ins.split(';'): arcpy.AddMessage(fc) arcpy.Clip_analysis(fc, clipfeats, out) Though I'm not entirely sure of ...


20

The hint provided by @mdsummer of using byid=TRUE works accurately. See the reproducible example, below: library(rgeos) library(sp) #Create SpatialPlygons objects polygon1 <- readWKT("POLYGON((-190 -50, -200 -10, -110 20, -190 -50))") #polygon 1 polygon2 <- readWKT("POLYGON((-180 -20, -140 55, 10 0, -140 -60, -180 -20))") #polygon 2 par(mfrow = c(1,...


19

After trying around with everything I finally figured out how to solve the problem. It had indeed to do with the CRS. Right click "Set CRS" was not enough here. I had to perform (on the raster) Raster->Projections->Warp, then set the desired CRS again and save as Geotiff. The mask layer (vector layer) had to be saved again with the same CRS. After that the ...


19

The main difference will be in the attributes of the results. When using Clip only the input feature’s attributes will be in the output (none from the clip feature), where if you used Intersect the attributes form all features used will be in the output.


19

The following quoted text came from a deleted blog post at infogeoblog.wordpress.com which was called Geo-Processing in QGIS. Given two input shapes: Clip creates a new shape based on the area of the input layer that is overlapped by the clipping layer. It is similar to the intersection but differs in that the attributes of the chosen layer only are ...


19

Another option here would be to perform a union and then delete the inner feature. After the union, select features from the initial layer and use the delete feature tool to remove the selected polygons.


19

I ran a test to determine how the speed and quality differs between the two methods, here are the results: Input data 4-band NAIP DOQQ image in .img format (349.34MB) A feature class used as the mask/clipper Performance Three trials were performed and benchmarked. The Clip (Data Management) method is significantly faster than the Extract by Mask (...


18

Dont use env.extent you need to get raster extent. import arcpy elevRaster = arcpy.sa.Raster('C:/data/elevation') myExtent = elevRaster.extent print myExtent i hope it helps you...


17

The following script clips polygon watersheds to polygon county boundaries, naming each output featureclass something like HspWBD_HU12_county name. Tested and it works. Make sure your values in the NAME field have no special characters or spaces (simple Python string methods can clean that up for you). import arcpy arcpy.env.workspace = r'D:\Projects\GDBs\...


17

It's not a very well-known feature, but you need the Feature Type Connections window. You access it like this (View > Windows > Feature Type Connections): In there select all the source feature types, select the transformer point to connect to, then click Connect:


17

st_intersection is probably the best way. Find whatever way works best to get an sf object to intersect with your input. Here's a way using the convenience of raster::extent and a mix of old and new. nc is created by example(st_read): st_intersection(nc, st_set_crs(st_as_sf(as(raster::extent(-82, -80, 35, 36), "SpatialPolygons")), st_crs(nc))) I don't ...


16

Since Erase (as @Jens linked) only is available with an Advanced license, you can download ET Geowizards. It can be installed as an Arcmap toolbox. Although you have to pay for the full suite, there's a free part of the program and the Erase function is included there (Overlay group).


15

you can use ogr2ogr if you want: ogr2ogr -clipsrc polygonforclipping.shp out.shp in.shp * * * -clipsrc [xmin ymin xmax ymax]|WKT|datasource|spat_extent: (starting with GDAL 1.7.0) clip geometries to the specified bounding box (expressed in source SRS), WKT geometry (POLYGON or MULTIPOLYGON), from a datasource or to the spatial extent of the -...


14

This question is similar to: Clip raster by raster with data extraction and resolution change, but coming from a different angle. However, I think the answer is likely the same. First off, choose which raster you wish to be definitive. I'll repeat my previous answer here for ease: Load required libraries: library(raster) library(rgdal) Read rasters: r1 = ...


14

In the Processing Toolbox, you can use the Clip vectors by extent tool from GDAL/OGR: Processing Toolbox > GDAL/OGR > [OGR] Geoprocessing > Clip vectors by extent


13

The "segments to approximate" option is what you are looking for. So the number you put into that filed will be the number of sides/quarter. So the default is 5, so you end up with a 20 sided polygon. Put in 25, and you end up with a 100 sided on, so it becomes smoother.


13

Since GDAL 2.1 (more info here) GDAL and OGR utilities can be used as library functions, so this task is incredibly simple now: from osgeo import gdal ds = gdal.Open('original.tif') ds = gdal.Translate('new.tif', ds, projWin = [-75.3, 5.5, -73.5, 3.7]) ds = None


12

If you are interested using Python, a good documentation is available at GeospatialPython.com, here. and clipraster.py source is here. The Process: Clipping a raster is a series of simple button clicks in high-end geospatial software packages. In terms of computing, geospatial images are actually very large, multi-dimensional arrays. Remote ...


12

You can use the dataset extent as a polygon geometry with the clip tool, as in the Using geometries in geoprocessing tools example. import arcpy pnt_array = arcpy.Array() extent = arcpy.Raster(in_raster).extent pnt_array.add(extent.lowerLeft) pnt_array.add(extent.lowerRight) pnt_array.add(extent.upperRight) pnt_array.add(extent.upperLeft) poly = arcpy....


12

ModelBuilder functions differently than batch processing in ArcGIS. Typically, you use iterators to loop through individual files rather than a spreadsheet-type list of files and actions, as in batch mode. The following is an example of the type of model you would need to loop through a workspace containing rasters in order to clip them to study area bound....


12

Seems to be a simple application of gDifference from the rgeos package: > require(rgeos) > ukhole = gDifference(uk, lnd) Warning message: In RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_difference") : spgeom1 and spgeom2 have different proj4 strings > plot(ukhole) The projection warning is because the LondonBoroughs shapefile doesn't have ...


11

If you want a solution that does not involve any extra extensions or "high-grade" licences you can try this: Union A1 and A2 to make A3. Then select by location where A3 does not have its center in A2 (you may need to select A3 where it DOES have its centre in A2 and then switch selection). OR select by attributes as Union will append attributes from A1 ...


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