Mapping

...

|

Maybe I’m a bad person for saying this, but I generally just accept that there is 1610 meters in a mile and 2.6 million square meters in a square mile. #roundingworks #gosueme

Math like this will probably be a bad thing if you are aiming a missile to hit a target thousands of miles away but for ordinary life, i think it’s quite fine. I use 1610 meters a mile, because it’s easier to round up — a half mile then becomes 805 meters, a quarter mile 402.5 miles, a tenth of a mile is 161 meters. If you rounded down to 1609 meters, that math becomes a lot harder to do in your head.
 
 
How did I not know that QGIS can natively open Shapefiles in ZIP files?
 
This is the best thing since like forever.

 

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.

Median Household Income in New York State, a Hexagram

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