<jupyter_start><jupyter_text>OGR Simple Features Library ogrmerge.py example<jupyter_code>%run "C:\Users\*USERNAME*\Anaconda3\pkgs\gdal-2.2.2-py36_1\Scripts\ogrmerge.py" -f GPKG -o merged.gpkg "C:\data\gdal\NE\10m_cultural\*.shp"<jupyter_output><empty_output><jupyter_text>OGR code examples<jupyter_code># Example 1 – create polygon geometry with OGR
from osgeo import ogr
r = ogr.Geometry(ogr.wkbLinearRing)
r.AddPoint(1,1)
r.AddPoint(5,1)
r.AddPoint(5,5)
r.AddPoint(1,5)
r.AddPoint(1,1)
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(r)
print(poly.ExportToWkt())
# Example 2 – create polygon geometry from GeoJSON
from osgeo import ogr
geojson = """{"type":"Polygon","coordinates":[[[1,1],[5,1],[5,5],[1,5], [1,1]]]}"""
polygon = ogr.CreateGeometryFromJson(geojson)
print(polygon) 
# Example 3 - basic geometric operations
# 1 create area
print("The area of our polygon is %d" % polygon.Area())
# 2 calculate centroid of polygon
cen = polygon.Centroid()
print(cen)
# 3 Get the boundary
b = polygon.GetBoundary()
print(b)
# 4 convex hull does the same in this case as boundary, as our polyogn is a square:
ch = polygon.ConvexHull() 
print(ch)
#5 buffer. A buffer value of 0 (zero) returns the same values as boundary and convex hull in this example:
buffer = polygon.Buffer(0) 
print(buffer)
# 6 check if a point is inside our polygon
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(10, 10)
polygon.Contains(point)
#Example 4 – writing polygon data to a newly created shapefile

import osgeo.ogr, osgeo.osr
# 1 set the spatial reference
spatialReference = osgeo.osr.SpatialReference()
spatialReference.ImportFromProj4('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')

# 2 create a new shapefile
driver = osgeo.ogr.GetDriverByName('ESRI Shapefile')
shapeData = driver.CreateDataSource('my_polygon.shp')

# 3 create the layer
layer = shapeData.CreateLayer('polygon_layer', spatialReference, osgeo.ogr.wkbPolygon)
layerDefinition = layer.GetLayerDefn()

# 4 geometry is put inside feature
featureIndex = 0
feature = osgeo.ogr.Feature(layerDefinition)
feature.SetGeometry(polygon)
feature.SetFID(featureIndex)

# 5 feature is put into layer
layer.CreateFeature(feature)
!ogrinfo my_polygon.shp
#Example 5 – using a spatial filter to select features

# import the modules
from osgeo import ogr
import os

# reference the shapefile and specify driver type
shapefile = r"C:\data\gdal\NE\10m_cultural\ne_10m_populated_places.shp"
driver = ogr.GetDriverByName("ESRI Shapefile")

# open the data source with driver, zero means open in read-only mode
dataSource = driver.Open(shapefile, 0)

# use the GetLayer() function for referencing the layer that holds the data
layer = dataSource.GetLayer()

# pass in the coordinates for the data frame to the SetSpatialFilterRect() function.
# this filter creates a rectangular extent and selects the features inside the extent
layer.SetSpatialFilterRect(-102, 26, -94, 36)
for feature in layer:
    # select only the cities inside of the USA
    # we can do this through a SQL query:
    # we skip the cities that are not in the USA,
    # and print the names of the cities that are
    if feature.GetField("ADM0NAME") != "United States of America":
        continue
    else:
        print(feature.GetField("NAME"))<jupyter_output><empty_output><jupyter_text>Shapely and Fiona<jupyter_code># Example 1 – create geometries with Shapely

from shapely.geometry import Polygon
p1 = Polygon(((1, 2), (5, 3), (5, 7), (1, 9), (1, 2)))
p2 = Polygon(((6,6), (7,6), (10,4), (11,8), (6,6)))
p1 
# A new command line is required for printing the second polygon:
p2
# Point takes tuples as well as positional coordinate values
from shapely.geometry import Point
point = Point(2.0, 2.0)
q = Point((2.0, 2.0))
q
# line geometry
from shapely.geometry import LineString
line = LineString([(0, 0), (10,10)])
line
# linear rings
from shapely.geometry.polygon import LinearRing
ring = LinearRing([(0,0), (3,3), (3,0)])
ring
# collection of points
from shapely.geometry import MultiPoint
points = MultiPoint([(0.0, 0.0), (3.0, 3.0)])
points
# collection of lines
from shapely.geometry import MultiLineString
coords = [((0, 0), (1, 1)), ((-1, 0), (1, 0))]
coords
# collection of polygons
from shapely.geometry import MultiPolygon
polygons = MultiPolygon([p1, p2,])
polygons
#Example 2 – apply geometrical methods with Shapely

print(p1.area)
print(p1.bounds)
print(p1.length)
print(p1.geom_type)
#Example 3 – reading JSON geometries with Shapely

import json
from shapely.geometry import mapping, shape
p = shape(json.loads('{"type": "Polygon", "coordinates": [[[1,1], [1,3 ], [3,3]]]}'))
print(json.dumps(mapping(p)))
p.area
# Example 4 – reading data with Fiona

import fiona
c = fiona.open(r"C:\data\gdal\NE\110m_cultural\ne_110m_admin_1_states_provinces.shp")
rec = next(iter(c))
rec.keys()
import pprint
pprint.pprint(rec['type'])
pprint.pprint(rec['id'])
pprint.pprint(rec['properties'])
pprint.pprint(rec['geometry'])
print(len(c))        # prints total amount of features     
print(c.driver)      # prints driver name
print(c.crs)         # prints coordinate reference system of data file
#Example 5 – accessing vector geometry in shapefiles using Shapely and Fiona

import pprint, fiona
with fiona.open\(r"C:\data\gdal\NE\110m_cultural\ne_110m_admin_1_states_provinces.shp") as src:
    pprint.pprint(src[0])
from shapely.geometry import shape
minnesota = {'type': 'Polygon', 'coordinates': [[(-89.61369767938538, 47.81925202085796), (-89.72800594761503, 47.641976019880644), (-89.84283098016755, 47.464725857119504), (-89.95765601272012, 47.286907253603175), (-90.13175391311144, 47.29274669045216), (-90.30585181350276, 47.29801768654593), (-90.47994971389409, 47.30385712339489), (-90.6540476142854, 47.309128119488676), (-90.85778194859611, 47.21282908791278), (-91.06097368036777, 47.117046820659795), (-91.26470801467849, 47.02126455340681), (-91.46844234898919, 46.9249655218309), (-91.59225908076053, 46.876260484395914), (-91.71661841507091, 46.82760712339301), (-91.8409777493813, 46.778385321635), (-91.96479448115262, 46.72970612241605), (-92.01189754918667, 46.71172272397848), (-92.27487891312, 46.65614472104858), (-92.26482784703924, 46.09522288673644), (-92.2965830146826, 46.09628225359842), (-92.54369971390233, 45.98569468849381), (-92.75696834997089, 45.88991242124082), (-92.89982784704182, 45.705763454768714), (-92.68922054723626, 45.51843638771068), (-92.76541744665064, 45.26708222104298), (-92.76647681351261, 44.9961426865367), (-92.79665584997105, 44.77602692318949), (-92.50507158076417, 44.58391978614445), (-92.38549231644078, 44.57492808692567), (-92.06215287959051, 44.43258535417769), (-91.949989183301, 44.364837551443046), (-91.87960588251947, 44.25742808692439), (-91.62770911331273, 44.085448920257036), (-91.28959021682704, 43.93729258887629), (-91.25729244664458, 43.854739488289965), (-91.2546569485977, 43.61397899024206), (-91.22819861526426, 43.5012468531974), (-92.54000484899348, 43.51977285417405), (-94.00102678096546, 43.51341665300214), (-95.35993608272871, 43.50018748633542), (-96.4526600817565, 43.50178945573647), (-96.43943091508982, 44.43576345476367), (-96.56061214881424, 45.39301768653837), (-96.73576941606754, 45.47081655535372), (-96.83470394569034, 45.62532908790642), (-96.78072791216147, 45.76079885515958), (-96.55689144568922, 45.87244578712617), (-96.53945064979075, 46.017966620460186), (-96.53890804725165, 46.199480088885366), (-96.60135901567638, 46.3513571233911), (-96.68548824744755, 46.51328522397769), (-96.73365068234358, 46.71647695574927), (-96.7458204821483, 46.9445250514533), (-96.77969438351563, 46.999043687521294), (-96.82041541216165, 47.29220408791309), (-96.82465287960957, 47.42661448830425), (-96.84423824744819, 47.54619375262766), (-96.89397681352912, 47.74886872007637), (-97.01515804725355, 47.954205023788134), (-97.13104244666809, 48.13729462339825), (-97.14850908078273, 48.31878225360731), (-97.16122148312654, 48.51458425556122), (-97.12734758175921, 48.64212169045756), (-97.1204746162644, 48.75852285419501), (-97.21413814979343, 48.90244171812793), (-97.22894344764504, 49.00088532164389), (-95.15883724646483, 48.9998259547819), (-95.15620174841791, 49.38401439065592), (-94.81754024939323, 49.38928538674975), (-94.64026424841592, 48.84001658791925), (-94.32912044958647, 48.670672919298795), (-93.63061011429733, 48.60928131773602), (-92.60984554723586, 48.45001455341253), (-91.63987891311746, 48.13993012144516), (-90.83026424840062, 48.27010305438836), (-89.59995174839577, 48.01027395282483), (-89.59940914585667, 48.01027395282483), (-89.4900319346622, 48.01340602289096), (-89.52269548211933, 47.96053538674391), (-89.61369767938538, 47.81925202085796)]]}
geom = shape(minnesota)
geom<jupyter_output><empty_output><jupyter_text>GeoPandas<jupyter_code># Example 1 – selecting and plotting geometry data with GeoPandas and Matplotlib

# import module
import geopandas as gpd
# magic command to use matplotlib plots inside Jupyter Notebook app
%matplotlib inline
# create a GeoDataFrame from shapefile and display it
df = gpd.read_file(r"C:\data\gdal\NE\110m_cultural\ne_110m_admin_1_states_provinces.shp")
df
# return object type
type(df)
# shape prints rows and columns
df.shape
# columns returns column names as a list item
df.columns
# subset data using pandas' dataframe .loc method
df.loc[0]
# print a list of all state names
df['name']
# reference a separate row and print all columns and attribute data
california = df.loc[df['name'] == "California"]
california
# plot state geometry
california.plot(figsize=(7,7))
# subset and plot multiple items using .iloc method
multipl = df.iloc[[5,7,9,11]]
multipl.plot(cmap="Set1", figsize=(7,7))
# subset items using .cx method and a bounding box and plot them
exp = df.cx[-124:-118,30:50]
exp.plot(cmap="Set1", figsize=(7,7))
# Example 2 – mapping wildfire data with GeoPandas
# import module
import geopandas
# import shapefile with state geometries
states = geopandas.read_file\(r"C:\data\gdal\NE\110m_cultural\ne_110m_admin_1_states_provinces.shp")
# display attribute table
states
# magic command for using matplotlib, next plot state geometries
%matplotlib inline
states.plot(figsize=(10,10))
# reference and display attribute table of shapefile with wildire data
fires = geopandas.read_file(r"C:\data\mtbs_fod_pts_data\mtbs_fod_pts_20170501.shp") 
fires
# plot all fires as point data on a map
fires.plot(markersize=1, figsize=(17,17))
# print coordinate reference system of firedata
fires.crs
# print coordinate reference system of state geometry shapefile
states.crs
# reproject fires shapefile so it´s equal to the state geometry shapefile
fires = fires.to_crs({'init': 'epsg:4326'})
# perform spatial join
state_fires = geopandas.sjoin(fires, states[['name', 'geometry']].copy(), op='within')
state_fires
# create pandas dataframe object with states and fire count
counts_per_state = state_fires.groupby('name').size()    
# list highest values first 
counts_per_state.sort_values(axis=0, ascending=False)
# add values to the original shapefile data as a new field
states = states.merge(counts_per_state.reset_index(name='number_of_fires'))        
# list first five rows of shapefile
states.head()
# create map with wildfire date per state
ax = states.plot(column='number_of_fires', figsize=(15, 6), cmap='OrRd', legend=True)
# same map with different colour scheme
ax = states.plot(column='number_of_fires', figsize=(15, 6), cmap='Accent', legend=True)
# same map, but without x and y axes and added title
import matplotlib.pyplot as plt
f, ax = plt.subplots(1, figsize=(18,6))
ax = states.plot(column='number_of_fires', cmap='Accent', legend=True, ax=ax)
lims = plt.axis('equal')
f.suptitle('US Wildfire count per state in 1984-2015')     
ax.set_axis_off()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13105726872246695, 'fraction_numerical': 0.23720772619451033, 'mean_word_length': 5.9196951934349356, 'pattern_counts': {'":': 5, '<': 17, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '333581', 'license': 'MIT', 'n_tokens_mistral': 6131, 'n_tokens_neox': 4388, 'n_words': 1055, 'path': 'Chapter05/Chapter5 Example Code.ipynb', 'repo_name': 'monocilindro/Mastering-Geospatial-Analysis-with-Python', 'size': '11804'}