Clip Vector Data to Area of Interest

Clip Geographic Shapes to AOI

Clip vector data in geographic coordinates to an AOI by using the geoclip function.

Read vector data from a shapefile of world cities into a geospatial table. The table represents the cities using point shapes in geographic coordinates. Extract the point shapes.

cities = readgeotable("worldcities.shp");
cShape = cities.Shape

cShape = 
  318×1 geopointshape array with properties:

               NumPoints: [318×1 double]
                Latitude: [318×1 double]
               Longitude: [318×1 double]
                Geometry: "point"
    CoordinateSystemType: "geographic"
           GeographicCRS: [1×1 geocrs]

Define a circular AOI with a radius of 11 degrees that is centered on Denver.

denver = geocode("Denver","city");
dAOI = aoicircle(denver,11);

Verify that the city data and the AOI use the same type of coordinate system.

isequal(cShape.CoordinateSystemType,dAOI.CoordinateSystemType)

ans = logical
   1

Verify that the city data and the AOI use the same geographic CRS.

isequal(cShape.GeographicCRS,dAOI.GeographicCRS)

ans = logical
   1

Clip the city data to the AOI.

gClippedG = geoclip(cShape,dAOI);

Display the AOI, the original city data, and the clipped city data on a map. Prepare to add a legend by specifying a display name for each plot.

figure
geobasemap grayland
geoplot(dAOI,FaceAlpha=0.2,DisplayName="AOI")
hold on
geoplot(cShape,"k",DisplayName="Original Data")
geoplot(gClippedG,"*",MarkerEdgeColor="#A2142F",DisplayName="Clipped Data")

Add a legend. Get a clearer view of the clipped data by zooming in.

legend
geolimits([18 56],[-130 -66])

Clip Projected Shapes to AOI

Clip vector data in projected coordinates to an AOI by using the mapclip function.

Read vector data from a shapefile of Boston roads into a geospatial table. The table represents the roads using line shapes in projected coordinates. Extract the line shapes.

roads = readgeotable("boston_roads.shp");
rShape = roads.Shape

rShape=2795×1 maplineshape array with properties:
                NumParts: [2795×1 double]
                Geometry: "line"
    CoordinateSystemType: "planar"
            ProjectedCRS: [1×1 projcrs]
      ⋮

Read a shapefile containing the coordinates of locations in Boston as a geospatial table. The table represents the locations using point shapes in projected coordinates. Create a subtable by geocoding four of the locations. Then, define a quadrangular AOI that bounds the locations.

locations = readgeotable("boston_placenames.shp");
sublocations = geocode(["BOSTON","BOSTON NECK","BEACON HILL","COPPS HILL"],locations);

bAOI = aoiquad(sublocations)

bAOI = 
  mappolyshape with properties:

              NumRegions: 1
                NumHoles: 0
                Geometry: "polygon"
    CoordinateSystemType: "planar"
            ProjectedCRS: [1×1 projcrs]

Verify that the road data and the AOI use the same type of coordinate system.

isequal(rShape.CoordinateSystemType,bAOI.CoordinateSystemType)

ans = logical
   1

Verify that the road data and the AOI use the same projected CRS.

isequal(rShape.ProjectedCRS,bAOI.ProjectedCRS)

ans = logical
   1

Clip the road data to the AOI.

rClipped = mapclip(rShape,bAOI);

Display the clipped road data and the boundary of the AOI on a map.

figure
geobasemap streets
geoplot(rClipped,LineWidth=1)
hold on
geoplot(bAOI,FaceColor="none",LineWidth=2)

Provide more geographic context for the data by zooming out.

geolimits([42.3494 42.3685],[-71.0774 -71.0461])