There are 2 basic types of geospatial data:

Vector data

Raster data

For vector data, there are 3 basic data types:

a. Point

b. LineString

c. Polygon

In [ ]:

```
from shapely.geometry import Point
```

In [ ]:

```
p1 = Point(13, 27.21)
p1
```

Out[ ]:

The code is importing the Point class from the shapely.geometry module. The Point class is a geometric object that represents a point in 2D space with x and y coordinates.

Then, the code creates a new Point object named p1 by calling the Point constructor with the arguments (13, 27.21). This creates a point at the coordinates (13, 27.21) in 2D space.

Similarly, we can create a second or multiple points.

In [ ]:

```
p2 = Point(61, 56.21)
p2
```

Out[ ]:

In [ ]:

```
p3 = Point(30, -6.21)
p3
```

Out[ ]:

In [ ]:

```
p4 = Point(2, -30.21)
p4
```

Out[ ]:

In [ ]:

```
from shapely.geometry import LineString
```

In [ ]:

```
line1 = LineString([p1, p2])
line1
```

Out[ ]:

The code creates a LineString object named line1 using the LineString class from the shapely.geometry module.

The LineString object represents a line in 2D space as a sequence of points. The LineString constructor takes a list of Point objects as its argument, so the code is passing a list containing two Point objects p1 and p2 as the argument to create a line segment connecting these two points.

In [ ]:

```
line2 = LineString([p1, p2, p3])
line2
```

Out[ ]:

In [ ]:

```
from shapely.geometry import Polygon
```

In [ ]:

```
Polygon([p1, p2, p3])
```

Out[ ]:

In [ ]:

```
Polygon([p1, p2, p3, p4])
```

Out[ ]:

We can have as many points as we need or want for polygons.

The `MultiLineString`

class in Python's Shapely library is used to represent collections of line objects in two-dimensional space. A MultiLineString object is a collection of one or more LineString objects that are considered as a single entity.

**Example:** dataset of the locations of weather stations in Bangladesh.

In geospatial data analysis, the **MultiPoint** class in Python's **Shapely** library is used to represent collections of point objects in two-dimensional space. A **MultiPoint** object is a collection of one or more **Point** objects that are considered as a single entity.

In [ ]:

```
from shapely.geometry import MultiPoint, MultiLineString, MultiPolygon
```

In [ ]:

```
multip1 = MultiPoint([p1, p2, p3, p4])
multip1
```

Out[ ]:

The `MultiLineString`

class in Python's Shapely library is used to represent collections of line objects in two-dimensional space. A `MultiLineString`

object is a collection of one or more `LineString`

objects that are considered as a single entity.

**Example:** dataset of roads in Dhaka.

The **MultiPolygon** class in Python's Shapely library is used to represent collections of polygon objects in two-dimensional space. A **MultiPolygon** object is a collection of one or more **Polygon** objects that are considered as a single entity.

**Example:** dataset of buildings footprints in New York City.

Let's create two polygons `poly1`

and `poly2`

from the four points we defined above and consider these polygons as a single entity -- a **MultiPolygon**.

In [ ]:

```
# imports the MultiPolygon class from the shapely.geometry module
from shapely.geometry import MultiPolygon
```

In [ ]:

```
poly1 = Polygon([p1, p2, p3])
poly2= Polygon([p2, p3, p4])
multi_poly = MultiPolygon([poly1, poly2])
multi_poly
```

Out[ ]:

In [ ]:

```
poly1
```

Out[ ]:

In [ ]:

```
poly2
```

Out[ ]:

In [ ]:

```
type(multi_poly)
```

Out[ ]:

With all the above objects, we can use the function `type`

to see or check their data type.

In [ ]:

```
p1.coords
```

Out[ ]:

In [ ]:

```
coordinates = list(p1.coords)
coordinates
```

Out[ ]:

We can see that the coordinates are given as a list of tuple. The first value is the x-cordinate and the second value is the y-cordiante. **Latitude is the x-cordinate and longitude is the y-coordinate.**

- First way using
`coords`

In [ ]:

```
# latitude
coordinates = list(p1.coords)
coordinates[0][0]
```

Out[ ]:

In [ ]:

```
# longitude
coordinates[0][1]
```

Out[ ]:

- Second way using
`x`

and`y`

attribute of the`Point`

object

In [ ]:

```
# latitude
p1.x
```

Out[ ]:

In [ ]:

```
# longitude
p1.y
```

Out[ ]:

- Using
`coords`

In [ ]:

```
line2.coords
```

Out[ ]:

In [ ]:

```
coordinates = list(line2.coords)
coordinates
```

Out[ ]:

It shows the coordinates of the points (here, three points) used to create this line.

To get the latitude of the third point, we can access it using the `coordinates`

list in the following way:

In [ ]:

```
coordinates[2][0]
```

Out[ ]:

- Using
`xy`

The `xy`

attribute of a `LineString`

object in Python's Shapely library returns a tuple containing two arrays representing the x and y coordinates of the line vertices.

In [ ]:

```
line2.xy
```

Out[ ]:

In [ ]:

```
# Getting longitude of the third point
line2.xy[1][2]
```

Out[ ]:

In [ ]:

```
poly1.exterior
```

Out[ ]:

In [ ]:

```
list(poly1.exterior.coords)
```

Out[ ]:

The `bounds`

attribute of a `Polygon`

or `MultiPolygon`

object in Python's Shapely library returns a tuple containing the minimum and maximum values of the x and y coordinates that define the bounding box of the polygon(s).

In [ ]:

```
line1.length
```

Out[ ]:

In [ ]:

```
poly1.length
```

Out[ ]:

In [ ]:

```
poly1.area
```

Out[ ]:

The centroid of `LineString`

or `Polygon`

is a point representing the "center of mass" or "average position" of the object. The centroid is a useful measure of the location and shape of the object.

In Shapely, the `centroid`

property of a `LineString`

or `Polygon`

object returns a `Point`

object that represents the centroid of the object. This `Point`

object has an x and y coordinate.

In [ ]:

```
line1.centroid
```

Out[ ]:

In [ ]:

```
print(line1.centroid)
```

In [ ]:

```
poly1.centroid
```

Out[ ]:

In [ ]:

```
print(poly1.centroid)
```

Distance from point 1 (`p1`

) to point 2 (`p2`

).

In [ ]:

```
p1.distance(p2)
```

Out[ ]:

Distance from point 3 (`p3`

) to point 1 (`p1`

).

In [ ]:

```
p3.distance(p1)
```

Out[ ]:

- Buffer around a
`Point`

object or geometry

In [ ]:

```
p1.buffer(33)
```

Out[ ]:

- Buffer around a
`MultiPoint`

object or geometry

In [ ]:

```
multip1.buffer(14)
```

Out[ ]:

- Buffer around a
`LineString`

object or geometry

In [ ]:

```
line1.buffer(12)
```

Out[ ]:

- Buffer around a
`Polygon`

object or geometry

In [ ]:

```
poly1.buffer(8)
```

Out[ ]:

In [ ]:

```
# just to see poly1 without buffer
poly1
```

Out[ ]:

The `contains`

method is a Boolean function from the Shapely library that returns `True`

if the point is inside the polygon, and `False`

otherwise.

In [ ]:

```
poly1.contains(p3) # p3 is one of the corners of this polygon
```

Out[ ]:

In [ ]:

```
new_point = Point(32, 28)
```

In [ ]:

```
# plot the Polygon and Point objects together
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.set_aspect('equal')
ax.plot(*poly1.exterior.xy)
ax.plot(new_point.x, new_point.y, 'o', color='red')
plt.show()
```

In [ ]:

```
poly1.contains(new_point)
```

Out[ ]:

Create a new point that is outside the polygon boundary.

In [ ]:

```
new_point2 = Point(25, 45)
```

In [ ]:

```
# plot the Polygon and Point objects together
fig, ax = plt.subplots()
ax.set_aspect('equal')
ax.plot(*poly1.exterior.xy)
ax.plot(new_point2.x, new_point2.y, 'o', color='red')
plt.show()
```

In [ ]:

```
poly1.contains(new_point2)
```

Out[ ]:

As expected, it says that this is out of the boundary of this polygon.

This is similar to `contains`

with the difference that it checks whether a object is within the boundary of another object.

In [ ]:

```
new_point.within(poly1)
```

Out[ ]:

In [ ]:

```
new_point2.within(poly1)
```

Out[ ]:

As the name suggests, `crosses`

allows us to see whether a geometry object (here, LineString) crosses another geometry object (here, a Polygon).

In [ ]:

```
from shapely.geometry import Point, Polygon, LineString
import matplotlib.pyplot as plt
# create multiple Point objects to be utilized for creating a Polygon object
p1 = Point(13, 27.21)
p2 = Point(61, 56.21)
p3 = Point(30, -6.21)
# use these points to create Polygon object
poly1 = Polygon([p1, p2, p3])
# create a Point object
point = Point(50, 6)
# create two LineString objects
line1 = LineString([(32, 28), (25, 45)])
line2 = LineString([(40, -5), (50, 20)])
# plot the Polygon, Point, and LineString objects together
fig, ax = plt.subplots()
ax.set_aspect('equal')
ax.plot(*poly1.exterior.xy)
ax.plot(point.x, point.y, 'o', color='red')
ax.plot(*line1.xy, color='blue')
ax.plot(*line2.xy, color='orange')
plt.show()
```

We see that `line1`

or blue line crosses the polygon. How can we check this without checking visually? The answer is `crosses`

.

In [ ]:

```
line1.crosses(poly1)
```

Out[ ]:

The result is as expected. Now, does the line `line2`

with orange color crosses the polygon? We know it doesn't, let's check this by using `crosses`

again.

In [ ]:

```
line2.crosses(poly1)
```

Out[ ]:

The `difference()`

method is used to compute the difference between two geometric objects. When this method is called on a Polygon object, it returns a new Polygon object that represents the difference between the original Polygon object and another Polygon object passed as an argument to the `difference()`

method.

In [ ]:

```
p1 = Point([20, 40])
```

In [ ]:

```
p2 = Point([40, 80])
```

In [ ]:

```
from shapely.geometry import Point
import matplotlib.pyplot as plt
# create two Point objects
point1 = Point(0.5, 0.5)
point2 = Point(1, 1)
# plot the Point objects together
fig, ax = plt.subplots()
ax.set_aspect('equal')
ax.plot(point1.x, point1.y, 'o', color='red')
ax.plot(point2.x, point2.y, 'o', color='blue')
plt.show()
```

In [ ]:

```
from shapely.geometry import Point, Polygon
import matplotlib.pyplot as plt
# create two Point objects
point1 = Point(0.5, 0.5)
point2 = Point(1, 1)
# create Polygon objects with buffer zone around each point
buffer1 = point1.buffer(0.6)
buffer2 = point2.buffer(0.6)
# plot the Polygon objects with buffer zone
fig, ax = plt.subplots()
ax.set_aspect('equal')
ax.fill(*buffer1.exterior.xy, alpha=0.5, color='red')
ax.fill(*buffer2.exterior.xy, alpha=0.5, color='blue')
plt.show()
```

In [ ]:

```
result = buffer1.difference(buffer2)
import matplotlib.pyplot as plt
from shapely.geometry import Point
# plot the difference between buffer1 and buffer2
fig, ax = plt.subplots()
ax.set_aspect('equal')
# ax.axis('off')
ax.fill(*result.exterior.xy, alpha=0.5, color='red')
plt.show()
```

In [ ]:

```
result2 = buffer2.difference(buffer1)
import matplotlib.pyplot as plt
from shapely.geometry import Point
# plot the difference between buffer2 and buffer1
fig, ax = plt.subplots()
ax.set_aspect('equal')
# ax.axis('off')
ax.fill(*result2.exterior.xy, alpha=0.5, color='blue')
plt.show()
```

In [ ]:

```
```

`intersection`

is used to obtain the region where two or more geometries overlap. When we apply the `intersection`

method to two geometries, the resulting geometry is a new geometry that represents the region where the two original geometries overlap. The intersection method can be applied to any type of geometry, including Points, Lines, and Polygons.

In [ ]:

```
intersection = buffer1.intersection(buffer2)
# plot the intersection between buffer1 and buffer2
fig, ax = plt.subplots()
ax.set_aspect('equal')
# ax.axis('off')
ax.fill(*intersection.exterior.xy, alpha=0.5, color='red')
plt.show()
```

We use `union`

to combine two or more geometries into a single geometry. When we apply the `union`

method to two geometries, it results into new geometry that encompasses the area covered by both of the original geometries. We can use the `union`

method with `Points`

, `Lines`

, and `Polygons`

.

In [ ]:

```
union = buffer1.union(buffer2)
# plot the union of buffer1 and buffer2
fig, ax = plt.subplots()
ax.set_aspect('equal')
# ax.axis('off')
ax.fill(*union.exterior.xy, alpha=0.5, color='red')
plt.show()
```

The `symmetric_difference()`

method in the Shapely library computes the symmetric difference between two geometric objects. It returns a new geometric object that represents the area of the geometries that do not overlap.

In [ ]:

```
union = buffer1.symmetric_difference(buffer2)
# plot the union of buffer1 and buffer2
fig, ax = plt.subplots()
ax.set_aspect('equal')
# ax.axis('off')
for polygon in union.geoms:
ax.fill(*polygon.boundary.xy, alpha=0.5, color='red')
plt.show()
```

We will be using GeoPandas which is basically a Pandas dataframe plus GeoSeries which is a column that contains vector geometry. These geometries are handled as Shapely object and thus the operations that we have already done for different Shapely object can also be applied to GeoPandas geometries.

First, download the shapefile from this link and upload it to your Google Drive. Dowload link: https://drive.google.com/drive/folders/1Z4XFVorvwiLVJ90Ph9fgbov2QHn4LJ09?usp=share_link

In [ ]:

```
from google.colab import drive
drive.mount('/content/drive')
```

In [ ]:

```
# Put the path to the shapefile in your Google Drive here
path = "/content/drive/MyDrive/Geospatial_Python/Data/"
```

In [ ]:

```
!pip install geopandas
```

In [ ]:

```
import geopandas as gpd
import matplotlib.pyplot as plt
# Read the shapefile
bd = gpd.read_file(path+'BGD_adm3.shp')
# Plot the shapefile
bd.plot()
plt.title("Upazilas in Bangladesh")
# Show the plot
plt.show()
```

We can also zoom in to the northeast part of Bangladesh by using `cx`

property of the GeoPandas dataframe. Here, we are specifying the range for latitude and longitude.

In [ ]:

```
bd.cx[91:92.5, 24:25.2].plot()
```

Out[ ]:

Let's have a look at the content of this GeoPandas dataframe.

In [ ]:

```
bd.head()
```

Out[ ]:

In the field of GIS, the columns are called attributes. Column `NAME_1`

represents the name of the divisions, column `NAME_2`

represents the name of the districts, and column `NAME_3`

represents the name of the upazilas.

Now, we can get the data (and map) of any division, district or upazila by using dataframe subsetting.

In [ ]:

```
# Get the data for Dhaka division and map it
dhaka = bd[bd['NAME_1'] == 'Dhaka']
```

In [ ]:

```
dhaka.head()
```

Out[ ]:

We can see that now we have data for Dhaka division only (note the values of `NAME_1`

column or attribute).

In [ ]:

```
dhaka.plot(color='green')
```

Out[ ]:

So, can we make a map of Rajshahi district? We can, by considering the field 'NAME_2' field corresponding to district name.

In [ ]:

```
rajshahi = bd[bd['NAME_2'] == 'Rajshahi']
rajshahi.plot()
```

Out[ ]:

Now, we will find out how many upazilas are in Sylhet division.

In [ ]:

```
sylhet_div = bd[bd['NAME_1'] == 'Sylhet'] # NAME_1 corresponds to division name
```

In [ ]:

```
sylhet_div.head()
```

Out[ ]:

In [ ]:

```
# Names of the upazilas
sylhet_div['NAME_3'].unique() # NAME_3 corresponds to upazila name
```

Out[ ]:

In [ ]:

```
# Number of upazila
len(sylhet_div['NAME_3'].unique())
```

Out[ ]:

In [ ]:

```
# Let's plot Bishwanath upazila
bishwanath = sylhet_div[sylhet_div['NAME_3'] == 'Bishwanath'] # subset the data to have data for Bishwanath upazila only
bishwanath.plot()
```

Out[ ]:

In [ ]:

```
# Let's plot Chhatak upazila under Sylhet division
sylhet_div[sylhet_div['NAME_3'] == 'Chhatak'].plot()
```

Out[ ]:

We can use `contextily`

package for adding basemaps.

In [ ]:

```
!pip install contextily
```

In [ ]:

```
import geopandas as gpd
import matplotlib.pyplot as plt
import contextily as ctx
# Set the figure size
fig, ax = plt.subplots(figsize=(10, 10))
# Plot the shapefile
bishwanath.plot(ax=ax, alpha=0.5)
# Add a basemap using contextily
ctx.add_basemap(ax=ax, source=ctx.providers.OpenStreetMap.Mapnik)
plt.title("The map of Bishwanath Upazila")
# Show the plot
plt.show()
```