Creating Album Covers with GIS—or Drawing Elevation as Vector Lines – xyHt
Mapping
Open Topograhpy – Infrastrucutre for Topographic Data
Open Topograhpy – Infrastrucutre for Topographic Data
9/17/21 by MapScaping
Web player: https://podcastaddict.com/episode/128610411
Episode: https://mcdn.podbean.com/mf/web/juwtpw/OpenTopography.mp3
OpenTopography provides community access to high-resolution, topography data, related tools, and resources … FOR FREE! OpenTopograhy.org Sponsored by Regrid.com Leading provider of land parcels and location context data for your maps, apps, and spatial analysis.
Convert Shapefiles into Hexagrams Automatically
Convert Shapefiles into Hexagrams Automatically
Mr. Bailey’s geogrid R-library is pretty neat. You will need to install the “R” language along with geogrid and sf libraries, which you can do with install.packages(‘geogrid’). After installing R and the library, the code to create a hexagram is pretty simple, although you will want to tweak the output a bit in QGIS after creating the computer-generated hexagram.
# required libraries
library(sf)
library(geogrid)
# load geojson shapefile
df <- st_read( "/tmp/nyco.geojson")
# iterate through seed values 1 to 10 to compare which
# is the best output of hex file (optional, but various seeds
# can improve the output)
for (i in seq(1,10)) {
# calculates hexgrid as a foreign object
tmp <- calculate_grid(shape = df, grid_type = "hexagonal", seed = i)
# Convert foreign object to an spatial dataframe
df_hex <- assign_polygons(df, tmp)
df_hex = st_as_sf(df_hex)
# save spatial dataframe as geojson
st_write(df_hex, paste("/tmp/nyco_hex_",i,".geojson", sep = "") )
}
Here is a sample hexagram I created using the Census TIGER/Line for New York counties.
How to create a spreadsheet converting or cross walking old Election Districts to new Election districts
How to create a spreadsheet converting or crosswalking old Election Districts to new Election districts 
It turns out with Python and GeoPANDAS available on pip it’s pretty easy to create a Microsoft Excel-style Spreadsheet with new election districts with the percentage of the old election districts within them. This code can also be used to calculate the percentage of a county in a Senate district. Obviously, you will need to update all of the paths and field names for your Shapefiles.
# requires the geopandas and pandas python libraries, available from pip
import pandas as pd
import geopandas as gpd
# Load Old EDs etc Shapefile
ac11 = gpd.read_file(r'/home/andy/Documents/GIS.Data/election.districts/2011albanyEDs.shp').to_crs(epsg='3857')
# Load the New EDs etc Shapefile
ac20 = gpd.read_file(r'/home/andy/Albany County October 2020 ED.gpkg').to_crs(epsg='3857')
# Save the area of new EDs
ac20['ed20_area']=ac20.area
# Create a GIS intersection of old and new EDs
edover=gpd.overlay(ac11, ac20, how='intersection')
# Calculate percentage of the New ED covered by the Old ED
edover['Percent of 2020 ED Area Covered by 2011 ED']=edover.area/edover['ed20_area']*100
# Rename columns to something more sensible
edover=edover.rename({'NAME_1':'2011 ED','NAME_2':'2020 ED'},axis=1)
# Exclude Portions of Old EDs that are less then 1% of New EDs
# Sort and Group by New ED and then Percentage of the Old ED
edover[(edover['Percent of 2020 ED Area Covered by 2011 ED']>1)][['2011 ED','2020 ED','Percent of 2020 ED Area Covered by 2011 ED']].sort_values(by=['2020 ED','Percent of 2020 ED Area Covered by 2011 ED']).groupby(by=['2020 ED','2011 ED']).median().to_excel('/tmp/New to Old Election District Crosswalk.xlsx')
Here is a sample of the Microsoft Excel-style spreadsheet that will be exported.
| 2020 ED | 2011 ED | Percent of 2020 ED Area Covered by 2011 ED |
| Albany Ward 1 Dist 1 | Albany 1-1 | 51.98 |
| Albany 8-4 | 47.97 | |
| Albany Ward 1 Dist 10 | Albany 2-7 | 100.00 |
| Albany Ward 1 Dist 2 | Albany 7-4 | 99.64 |
| Albany Ward 1 Dist 3 | Albany 1-1 | 53.42 |
| Albany 1-2 | 46.35 | |
| Albany Ward 1 Dist 4 | Albany 1-3 | 46.06 |
| Albany 1-4 | 53.77 | |
| Albany Ward 1 Dist 5 | Albany 1-1 | 100.00 |
| Albany Ward 1 Dist 6 | Albany 1-5 | 96.07 |
| Albany 2-2 | 3.72 | |
| Albany Ward 1 Dist 7 | Albany 1-1 | 28.87 |
| Albany 2-10 | 34.49 | |
| Albany 2-8 | 36.52 | |
| Albany Ward 1 Dist 8 | Albany 1-5 | 12.52 |
| Albany 1-6 | 30.09 | |
| Albany 2-2 | 24.88 | |
| Albany 2-7 | 32.51 | |
| Albany Ward 1 Dist 9 | Albany 1-4 | 92.60 |
| Albany 1-5 | 7.09 | |
| Albany Ward 10 Dist 1 | Albany 10-2 | 99.30 |
| Albany Ward 10 Dist 2 | Albany 10-8 | 99.69 |
| Albany Ward 10 Dist 3 | Albany 10-1 | 99.64 |
| Albany Ward 10 Dist 4 | Albany 10-4 | 99.57 |
| Albany Ward 10 Dist 5 | Albany 10-5 | 99.44 |
| Albany Ward 10 Dist 6 | Albany 10-5 | 68.17 |
| Albany 10-8 | 31.83 | |
| Albany Ward 10 Dist 7 | Albany 10-8 | 48.18 |
| Albany 11-5 | 13.05 | |
| Albany 11-7 | 38.66 | |
| Albany Ward 10 Dist 8 | Albany 11-7 | 32.47 |
| Albany 6-1 | 67.24 |
A Better Python GeoEncoding Script to Use with NYS SAM
I rewrote the short little geocoder that I posted last night in Python, using the multiple address function in NY’s Street and Address Maintenance (SAM) Program. This should be able to encode up to 1,000 lines in a text file. You would use this script by typing python sam-geocode.csv with one address per line.
By running python sam-geocode.py address.txt, with a file like this:
5 Winding Brook Drive Guilderland NY 12084
1224 Union Street Schenectady NY 12308
6239 Randomwood Drive Schenectady NY 12030
1529 Western Avenue Suite 102 Albany NY 12203
4017b State Street Schenectady NY 12304
4017b State Street Schenectady NY 12304
1529 Western Avenue Suite 102 Albany NY 12203
1529 Western Avenue Suite 102 Albany NY 12203
2842 W. Lydius Street Schenectady NY 12306
200 Bloomingdale Road Altamont NY 12009
Gives you out a CSV file named address-geocoded.csv.
"1224 Union St, Schenectady, NY, 12308",42.81103508600006,-73.92220575399995
"6239 Randomwood Dr, Schenectady, NY, 12303",42.74713565800005,-73.93299995099994
"1529 Western Ave, Albany, NY, 12203",42.68293244500006,-73.84222237599994
"4017 B State St, Schenectady, NY, 12304",42.76677680800003,-73.88794539599996
"4017 B State St, Schenectady, NY, 12304",42.76677680800003,-73.88794539599996
"1529 Western Ave, Albany, NY, 12203",42.68293244500006,-73.84222237599994
"1529 Western Ave, Albany, NY, 12203",42.68293244500006,-73.84222237599994
"2842 W Lydius St, Schenectady, NY, 12306",42.75003186400005,-73.95573429299998
"23 Dresden Ct, Albany, NY, 12203",42.69198706100008,-73.87999566599996
"629 Salvia Ln, Schenectady, NY, 12303",42.73869920900006,-73.93268336399996
"5 Winding Brook Dr, Guilderland, NY, 12084",42.69401880400005,-73.90855096599995
"3 Feeney Ln, Schenectady, NY, 12303",42.72653497400006,-73.90924777399994
"730 Sachem Cir, Slingerlands, NY, 12159",42.67422884300004,-73.88860130999996
"200 Bloomingdale Ln, Altamont, NY, 12009",42.69634668300006,-73.94903687499993
Now, I could improve error capture on this script but for my purposes this is very simple and works wells for well-formed addresses. Maybe not as simple as the PHP script I posted last night, but this way of Geocoding is much faster and written in a language that’s more appropriate for such purposes.
#!/usr/bin/python
import requests,sys,json,csv,os
query = '{"records": ['
f = open(sys.argv[-1], 'r+')
with open(sys.argv[-1]) as f:
i = 0
for line in f:
query += '{ "attributes": { "OBJECTID":'+str(i)+', "SINGLELINE": "'+line.rstrip()+'"} },'+"n"
i+=1
query += ']}'
post = { 'f':'pjson', 'outSR': 4326, 'addresses': query }
url = 'https://gisservices.its.ny.gov/arcgis/rest/services/Locators/Street_and_Address_Composite/GeocodeServer/geocodeAddresses'
req = requests.post(url, data = post)
addresses = json.loads(req.text)['locations']
csv = open(os.path.splitext(sys.argv[-1])[0]+'-geocode.csv', "w")
for place in addresses:
csv.write('"'+place['address'].replace('"', '\"')+'",'+str(place['location']['y'])+','+str(place['location']['x'])+"n")
csv.close()
This also works for most other states, as most states have geocoding endpoints, as otherwise the fire department wouldn’t necessarily be able to find your house when it’s on fire when you call 911. To find your state, try searching for StateName "GeocodeServer/geocodeAddresses". For example, replace URL with:
Massachusetts – https://gis.massdot.state.ma.us/arcgis/rest/services/General/MassDOTStreetAddressLocator/GeocodeServer/geocodeAddresses
Pennsylvania – https://maps.pasda.psu.edu/arcgis/rest/services/apps/AddressLocator/GeocodeServer/geocodeAddresses
Spatial join — Intro to Python GIS documentation
Geopandas
Geopandas! 🗺
Teaching myself about the Geopandas library. Not because I’m interested in the plotting capacities but it sure would be nice to have a fast “geocoder” that could join point data to other shapefiles I have like counties, municipalities, Assembly districts, etc. without using QGIS.
For example, I can feed the state’s address management geocoder up to 1,000 addresses and it will spit back a list of coordinates. In theory quite easily I could take those coordinates, spatially join them with Geopandas and add a column with the Assembly district and municipality – all without loading the shapefile or plotting a map or doing scripting within QGIS.
https://automating-gis-processes.github.io/CSC18/lessons/L4/spatial-join.html