# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# ## GeoPandas Demo: Get Counties
# This example demonstrates how to grab data from an ArcGIS MapService and pull it into a GeoPandas data frame.

# +
import requests
import pandas as pd
import geopandas as gpd

# %matplotlib inline
# -

# ### Fetching some data
# We'll tap into a NOAA map server to pull some state boundary features...
# * Build the request
# * Send the request, receive the response

# +
#Build the request and parameters to fetch county features
#  from the NOAA ArcGIS map server end point
stateFIPS = '37' #This is NC

url = 'https://nowcoast.noaa.gov/arcgis/rest/services/nowcoast/mapoverlays_political/MapServer/find'
params = {'searchText':stateFIPS,
          'contains':'true',
          'searchFields':'STATEFP',
          'sr':'',
          'layers':'2',
          'layerDefs':'',
          'returnGeometry':'true',
          'maxAllowableOffset':'',
          'geometryPrecision':'',
          'dynamicLayers':'',
          'returnZ':'false',
          'returnM':'false',
          'gdbVersion':'',
          'returnUnformattedValues':'false',
          'returnFieldName':'false',
          'datumTransformations':'',
          'layerParameterValues':'',
          'mapRangeValues':'',
          'layerRangeValues':'',
          'f':'json'}
# -

#Fetch the data
response = requests.get(url,params)

# ### Examining the response
# * Convert the response to a JSON object
# * Examine its structure
# * Extract the `attributes` and `geometry` elements.

#Convert to a JSON object (i.e. a dictionary)
respons_js = response.json()

#The 'results' object contains a record for each county returned, i.e., a feature
results = respons_js['results']
len(results)

#Within each item in the results object are the following items
results[0].keys()

#The 'attributes' item contains the feature attributes
results[0]['attributes']

#And the geometry object contains the shape
results[0]['geometry']

# ### Convert the elements to dataFrames
# * Creating a dataFrame from the Results object
# * "Exploding" the dictionary values in the `attributes` and `geometry` columns
# * Concatenating dataFrames lengthwise (adding columns)

#Create a dataFrame from the results, 
#  keeping just the attributes and geometry objects
df = pd.DataFrame(results,columns=('attributes','geometry'))
df.head()

#Explode the dictionary values into fields
dfCounties = df['attributes'].apply(pd.Series)
dfGeom = df['geometry'].apply(pd.Series)

#Combine the two
dfAll = pd.concat((dfCounties,dfGeom),axis='columns')
dfAll.head()

# ### Converting the [ESRI] geometry coordinates to a [shapely] geometric feature
# The `dfAll` dataframe now has all feature attributes and the geometry object stored in the `rings` column. 
# * Exploring the 'rings' object
# * Exploring the `shapely` package: rings, polygons, and multipolygons
# * Using shapely to create features
# * Converting the dataFrame to geodataFrame
# * Plotting the output

#Explore the values in the "ring" column, looking at the first row of data
rings = dfAll['rings'][0]
print ("There is/are {} ring(s) in the record".format(len(rings)))
print ("There are {} vertices in the first ring".format(len(rings[0])))
print ("The first vertex is at {}".format(rings[0][0]))

# So, the "ring" value in each row of our dataframe contains a *list* of rings, with each ring being a list of coordinates defining the vertices of our polyon. Usually the list of rings only includes one ring, the outer boundary of a single polygon. However, it's possible it contains more than one, e.g. the boundary of Hawaii. 

# Now we'll extract the first ring object from the ring list of the first record in our dataframe and convert it to a Shapely polygon object. To do this we need to import a few Shapely geometry class objects.

#Import the shapely objects we'll need
from shapely.geometry import LinearRing
from shapely.geometry import Polygon

#Create a shapely polygon from the first ring 
ring = rings[0]        # Get the outer ring, in coordinates
r = LinearRing(ring)   # Convert coordinates to shapely ring object
s = Polygon(r)         # Convert shapely ring object to shapely polygon object
s.area                 # Show the area of the polygon


# Now that we've seen the proof of concept, we'll form a Python function that 
# * takes a list of rings (i.e., the value of one row's `rings` field), 
# * converts each ring item in this ring list into a Shapely LinearRing object, 
# * converts *that* into a Shapely polygon object, adding each these polygons to a list, 
# * and then constructs a Shapely MultiPolygon object from the list of polygons

#A function to convert all rings into a Shapely multipolygon object
def polyFromRings(rings):
    #Import necessary Shapely classes
    from shapely.geometry import LinearRing, Polygon, MultiPolygon
    #Construct an empty list of polygons
    polyList = []
    #Compile a list of shapely ring objects and convert to polygons
    for ring in rings:
        #Construct a ring from the ring coordinates
        r = LinearRing(ring)
        #Convert the ring to a shapely polygon
        s = Polygon(r)
        #Add the polygon to the polyList
        polyList.append(s)
    #Convert the list of polyongs to a multipolygon object
    multiPoly = MultiPolygon(polyList)
    return multiPoly


# Now, we use Panda's `apply` method to apply the "polyFromRings" function above to each row's "ring" values. 

#Apply the function to each item in the geometry column
dfAll['geometry']=dfAll.apply(lambda x: polyFromRings(x.rings),axis='columns')

# ### Convert dataframe to a *geo*dataframe
# With the rings successfuly converted to Shapely geometry objects, we can now "upgrade" our Pandas dataframe to a GeoPandas dataframe, capable of spatial analysis

#Create a geodataframe from our pandas dataframe (the geometry column must exist)
gdf=gpd.GeoDataFrame(dfAll)
gdf.head()

#Set the projection (obtained from spatialReference column)
gdf.crs = {'init': 'epsg:3857'}

#Check the data types; note some should be fixed!
gdf.dtypes

#Convert the `ALAND` and `AWATER` to floating point values
gdf['ALAND']=gdf['ALAND'].astype('double')
gdf['AWATER']=gdf['AWATER'].astype('double')
gdf.dtypes

#Use familiar Pandas operation to select a feature and 
gdf[gdf['NAME'] == 'Durham'].plot();

# ### Saving the data
# We can save the attributes to CSV file or save the feature class to a shapefile

#Save our attribute data to a shapefile
gdf.to_csv("counties_{}.csv".format(stateFIPS))

# Saving data to a shapefile is a bit more finicky. In particular, we need to remove the old "rings" field in the geodataframe.

#Delete the 'rings' column from our geodataframe
gdf.drop('rings',axis='columns',inplace=True)

#Write the data to a file
gdf.to_file(driver='ESRI Shapefile',filename="./data/NC_Counties.shp")

# ## Recap
# So here we've imported a layer from an ESRI Map Service and done the necessary conversions to get this into a GeoPandas dataframe -- and also export it. This reveals a bit about the requirements of a geodataframe, namely the structure of the geometry column and how the Shapely package helps with that. 
#
# From here, we can explore more about the cool things we can do with a GeoPandas dataframe.

StarcoderdataJupyterScriptsDedupFiltered