What is Geospatial Data?¶
Overview¶
- Points - Combination of an X and Y coordinate. Longitude and Latitude in geospatial terms.
- Lines - One or more connections between Points.
- Polygons - Plane figure that is described by a finite number of straight line segments connected to form a closed polygonal chain. Typically a collection of Points.
- Multipolygons - Combination of Polygons.
- Rasters - Raster graphic represents a two-dimensional picture as a rectangular matrix or grid of square pixels.
Points¶
In [1]:
from shapely.geometry import Point
point = Point(0.0, 0.0)
print("Area", point.area)
print("Length", point.length)
print("Bounds", point.bounds)
print("Longitude", point.x)
print("Latitude", point.y)
Area 0.0 Length 0.0 Bounds (0.0, 0.0, 0.0, 0.0) Longitude 0.0 Latitude 0.0
In [2]:
import matplotlib.pyplot as plt
points = [Point(-4, 0.),
Point(-5, 12.),
Point(-6., 3.),
Point(-8., 7.),
Point(-4., 8.),
Point(-3., 2.),
Point(-1., 6.)]
xs = [point.x for point in points]
ys = [point.y for point in points]
plt.scatter(xs, ys)
plt.show()
LineStrings¶
In [3]:
from shapely.geometry import LineString
line = LineString([(0, 0), (1, 1)])
print("Area", line.area)
print("Length", line.length)
print("Bounds", line.bounds)
print("Longitude, Latitude", line.xy)
Area 0.0
Length 1.4142135623730951
Bounds (0.0, 0.0, 1.0, 1.0)
Longitude, Latitude (array('d', [0.0, 1.0]), array('d', [0.0, 1.0]))
In [4]:
Lines = [LineString([Point(-4, 0.), Point(-5, 12.), Point(-6., 3.)]),
LineString([Point(-8., 7.),Point(-4., 8.),Point(-3., 2.),Point(-1., 6.)])]
plt.plot(Lines[0].xy[0], Lines[0].xy[1])
plt.plot(Lines[1].xy[0], Lines[1].xy[1])
plt.show()
Polygons¶
In [5]:
from shapely.geometry import Polygon
polygon = Polygon([(0, 0), (1, 1), (1, 0)])
print("Area", polygon.area)
print("Length", polygon.length)
print("Bounds", polygon.bounds)
print("Boundary", polygon.boundary)
print("Center", polygon.centroid)
Area 0.5 Length 3.414213562373095 Bounds (0.0, 0.0, 1.0, 1.0) Boundary LINESTRING (0 0, 1 1, 1 0, 0 0) Center POINT (0.6666666666666666 0.3333333333333333)
In [6]:
Polygons = Polygon([(-4, 0.), (-5, 0.),(-6., 3.), (-8., 7.),(-4., 8.),(-3., 2.), (-1., 2.)])
plt.figure(figsize=(10, 5))
plt.fill(*Polygons.exterior.xy)
plt.show()
In [7]:
import shapely.ops as so
import matplotlib.pyplot as plt
from itertools import combinations
r1 = Polygon([(0,0),(0,1),(1,1),(1,0),(0,0)])
r2 = Polygon([(1,1),(1,2),(2,2),(2,1),(1,1)])
r3 = Polygon([(0.5,0.5),(0.5,1.5),(1.5,1.5),(1.5,0.5),(0.5,0.5)])
r4 = Polygon([(3,3),(3,4),(4,4),(4,3),(3,3)])
polygons = [r1, r2, r3, r4]
new_shape = so.unary_union([r1, r2, r3])
multipolygon = so.unary_union([r1, r2, r3, r4])
In [8]:
import geopandas as gpd
import pandas as pd
points = [Point(-4, 0.), Point(-5, 12.),Point(-6., 3.),Point(-8., 7.),
Point(-4., 8.),Point(-3., 2.),Point(-1., 6.)]
Lines = [LineString([Point(-4, 0.), Point(-5, 12.),
Point(-6., 3.)]),
LineString([Point(-8., 7.),Point(-4., 8.),
Point(-3., 2.),Point(-1., 6.)])]
geoms = Polygon([(0,0),(0,5),(5,5),(5,0),(0,0)])
In [9]:
points_gdf = gpd.GeoDataFrame(pd.DataFrame({'Name': ["P1","P2","P3","P4","P5","P6","P7"]}),
crs="EPSG:4326",
geometry=points)
lines_gdf = gpd.GeoDataFrame(pd.DataFrame({'Name': ["Line 1", "Line 2"]}),
crs="EPSG:4326",
geometry=Lines)
polygon_gdf = gpd.GeoDataFrame(pd.DataFrame({'Name': ["Polygon"]}),
crs="EPSG:4326",
geometry=[geoms])
gdf = pd.concat([points_gdf, lines_gdf, polygon_gdf]).reset_index(drop=True)
gdf["area"] = gdf.area
gdf['boundary'] = gdf.boundary
gdf['centroid'] = gdf.centroid
C:\Users\AdamS\AppData\Local\Temp\ipykernel_21376\892064037.py:14: UserWarning: Geometry is in a geographic CRS. Results from 'area' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation. gdf["area"] = gdf.area C:\Users\AdamS\AppData\Local\Temp\ipykernel_21376\892064037.py:16: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation. gdf['centroid'] = gdf.centroid
In [10]:
print(gdf[['Name', 'area', 'boundary', 'centroid']])
Name area boundary \
0 P1 0.0 GEOMETRYCOLLECTION EMPTY
1 P2 0.0 GEOMETRYCOLLECTION EMPTY
2 P3 0.0 GEOMETRYCOLLECTION EMPTY
3 P4 0.0 GEOMETRYCOLLECTION EMPTY
4 P5 0.0 GEOMETRYCOLLECTION EMPTY
5 P6 0.0 GEOMETRYCOLLECTION EMPTY
6 P7 0.0 GEOMETRYCOLLECTION EMPTY
7 Line 1 0.0 MULTIPOINT (-4.00000 0.00000, -6.00000 3.00000)
8 Line 2 0.0 MULTIPOINT (-8.00000 7.00000, -1.00000 6.00000)
9 Polygon 25.0 LINESTRING (0.00000 0.00000, 0.00000 5.00000, ...
centroid
0 POINT (-4.00000 0.00000)
1 POINT (-5.00000 12.00000)
2 POINT (-6.00000 3.00000)
3 POINT (-8.00000 7.00000)
4 POINT (-4.00000 8.00000)
5 POINT (-3.00000 2.00000)
6 POINT (-1.00000 6.00000)
7 POINT (-4.92923 6.64384)
8 POINT (-3.74523 5.39758)
9 POINT (2.50000 2.50000)
In [11]:
geoms = [Polygon([(0,0),(0,5),(5,5),(5,0),(0,0)]),
Polygon([(2.5,2.5),(2.7,7.5),(7.5,7.5),(7.5,2.5),(2.5,2.5)])]
polygon1_gdf = gpd.GeoDataFrame(pd.DataFrame({'Name': ["Polygon1"]}),
crs="EPSG:4326",
geometry=[geoms[0]])
polygon2_gdf = gpd.GeoDataFrame(pd.DataFrame({'Name': ["Polygon1"]}),
crs="EPSG:4326",
geometry=[geoms[1]])
In [12]:
intersection = polygon2_gdf.overlay(polygon1_gdf, how='intersection')
union = polygon2_gdf.overlay(polygon1_gdf, how='union')
symetrical_difference = polygon2_gdf.overlay(polygon1_gdf, how='symmetric_difference')
difference = polygon2_gdf.overlay(polygon1_gdf, how='difference')
In [13]:
gdf = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
gdf.plot()
plt.show()
Rasters¶
In [14]:
import rasterio
dataset = rasterio.open('../resources/assets/GC.jpg')
print("Width", dataset.width)
print("Height", dataset.height)
print("Bounds", dataset.bounds)
print("CRS", dataset.crs)
plt.imshow(dataset.read(1))
plt.show()
Width 1024 Height 1024 Bounds BoundingBox(left=0.0, bottom=1024.0, right=1024.0, top=0.0) CRS None
C:\Users\AdamS\anaconda3\envs\pythonmaps-tutorial\lib\site-packages\rasterio\__init__.py:220: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix be returned. s = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)
Projections¶
In [15]:
import matplotlib.ticker as mticker
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
plt.figure(figsize=(10, 5))
ax = plt.axes(projection=ccrs.Mercator())
gdf.plot(ax=ax, transform=ccrs.PlateCarree())
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=2, color='black', alpha=0.5, linestyle='--')
plt.show()
In [16]:
plt.figure(figsize=(10, 5))
ax = plt.axes(projection=ccrs.LambertCylindrical())
gdf.plot(ax=ax, transform=ccrs.PlateCarree())
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=2, color='black', alpha=0.5, linestyle='--')
plt.show()
