Query Points in Area of Interest

Since R2024b

This topic shows how to determine if points are within an area of interest (AOI). Within the topic, you:

• Specify the points to query using point shape objects.

• Define the AOI. Mapping Toolbox™ represents AOIs using polygon shape objects. For information about creating AOIs, see Define Areas of Interest.

• Determine whether the locations are within the AOI by using the isinterior function.

To query points within an AOI, these properties of the point shapes and the AOI must match:

• The type of coordinate system. You can get the type of coordinate system for a shape object by querying the CoordinateSystemType property.

• The coordinate reference system (CRS). For shape objects in a geographic coordinate system, you can get the CRS by querying the GeographicCRS property. For shape objects in a planar coordinate system, you can get the CRS by querying the ProjectedCRS property.

These examples show how to determine whether geographic and projected point objects are within an AOI.

Query Geographic Shapes in AOI

Geocode the names of several European cities. The result is a geospatial table that represents the cities using point shapes in geographic coordinates. Extract the point shapes.

names = ["Madrid","Dublin","London","Berlin","Rome", ...
    "Prague","Budapest","Milan","Amsterdam","Munich"];
cities = geocode(names,"city");
cShape = cities.Shape

cShape = 
  10×1 geopointshape array with properties:

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

Geocode the city of Paris. Define a circular AOI with a radius of 6.5 degrees that is centered on Paris.

paris = geocode("Paris","city");
parisAOI = aoicircle(paris,6.5)

parisAOI = 
  geopolyshape with properties:

              NumRegions: 1
                NumHoles: 0
                Geometry: "polygon"
    CoordinateSystemType: "geographic"
           GeographicCRS: [1×1 geocrs]

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

isequal(cShape.CoordinateSystemType,parisAOI.CoordinateSystemType)

ans = logical
   1

Verify that the cities and the AOI use the same geographic CRS.

isequal(cShape.GeographicCRS,parisAOI.GeographicCRS)

ans = logical
   1

Determine if the cities are within the AOI. The isinterior function returns an array of logical values, where 1 indicates that the corresponding city is within the AOI, and 0 indicates that the corresponding city is not within the AOI.

inParisAOI = isinterior(parisAOI,cShape)

inParisAOI = 10×1 logical array

   0
   0
   1
   0
   0
   0
   0
   1
   1
   1

Create two arrays. The first array contains the cities that are in the AOI. The second array contains the cities that are not in the AOI.

citiesInAOI = cShape(inParisAOI);
citiesOutAOI = cShape(~inParisAOI);

Display the AOI, the cities that are in the AOI, and the cities that are not in the AOI on a map. Prepare to add a legend by specifying a display name for each plot.

figure
geobasemap grayland
geoplot(parisAOI,FaceAlpha=0.2,DisplayName="AOI")
hold on
geoplot(citiesInAOI,"o",MarkerFaceColor="#77AC30",MarkerEdgeColor="k", ...
    DisplayName="Inside")
geoplot(citiesOutAOI,"x",MarkerEdgeColor="#A2142F",MarkerSize=8, ...
    DisplayName="Outside")

Add a legend. Provide more geographic context for the data by zooming out.

legend
geolimits([30 63],[-26 38])

Query Projected Shapes in AOI

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. Extract the point shapes.

locations = readgeotable("boston_placenames.shp");
shape = locations.Shape

shape = 
  13×1 mappointshape array with properties:

               NumPoints: [13×1 double]
                       X: [13×1 double]
                       Y: [13×1 double]
                Geometry: "point"
    CoordinateSystemType: "planar"
            ProjectedCRS: [1×1 projcrs]

Create a subtable by geocoding two of the locations. Then, create an AOI of circles from the locations. Specify the radius of each circle as 700 meters.

sublocations = geocode(["BEACON HILL","BROAD CANAL"],locations);
aoi = aoicircle(sublocations,700)

aoi = 
  mappolyshape with properties:

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

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

isequal(shape.CoordinateSystemType,aoi.CoordinateSystemType)

ans = logical
   1

Verify that the locations and the AOI use the same projected CRS.

isequal(shape.ProjectedCRS,aoi.ProjectedCRS)

ans = logical
   1

Determine if the locations are within the AOI. The isinterior function returns an array of logical values, where 1 indicates that the corresponding location is within the AOI, and 0 indicates that the corresponding location is not within the AOI.

inAOI = isinterior(aoi,shape)

inAOI = 13×1 logical array

   0
   0
   1
   1
   0
   1
   0
   0
   1
   1
      ⋮

Create two arrays. The first array contains the locations that are in the AOI. The second array contains the locations that are not in the AOI.

locInAOI = shape(inAOI);
locOutAOI = shape(~inAOI);

Display the AOI, the locations that are in the AOI, and the locations that are not in the AOI on a map. Prepare to add a legend by specifying a display name for each plot.

figure
geoplot(aoi,FaceAlpha=0.2,DisplayName="AOI")
hold on
geoplot(locInAOI,"o",MarkerFaceColor="#77AC30",MarkerEdgeColor="k", ...
    DisplayName="Inside")
geoplot(locOutAOI,"x",MarkerEdgeColor="#A2142F", ...
    DisplayName="Outside")

Add a legend.

legend