Another Example of GeoPandas and Python Spatial Joins πΊ
Yesterday, I posted a much more complicated piece of code that pulled addresses from the SAM (State Address Management) database and did a spatial join to add a column to the file with Assembly District and Municipality. This was a bit too complex, so I made a simpler one for other purposes that doesn’t require the coordinates to be obtained from SAM.
This python script takes two parameters:
Path to a CSV file that contains an X and Y coordinate
Path to a Shapefile or Geopackage to Join Against
Then the code will create a new CSV file with the spatially joined attributes pulled from the Shapefile. I have only run it on a few large data sets, but I found it took roughly 1 second to join 1,000 records from call to end of end of script.
#!/usr/bin/python
import requests,sys,json,os,csv
import pandas as pd
import geopandas as gpd
lines=[]
# read list of addresses from parameter 1
with open(sys.argv[-2], newline='') as csvfile:
for line in csv.DictReader(csvfile):
lines.append(line)
# convert to pandas
locPd = pd.DataFrame(lines,columns=lines[0].keys())
locPd.convert_dtypes()
locPd = gpd.GeoDataFrame(locPd, geometry=gpd.points_from_xy(locPd.x.astype('float32'), locPd.y.astype('float32')))
# run spatial joins against parameter 2
ad = gpd.read_file(sys.argv[-1])
locPd = gpd.sjoin(locPd, ad, op="within")
# remove added geometery and index columns
del locPd['geometry']
del locPd['index_right']
# write pandas back to out csv
locPd.to_csv (os.path.splitext(sys.argv[-2])[0]+'-output.csv', index = False, header=True)
It seems like most of the state’s GIS services are currently down.
Well, the ones at Building 5 in State Office Campus like SAM, orthphotography, tax maps, etc … the DEC server is still working, although that’s probably because they are hosted in ESRI cloud.
I was reading the state GIS’ email list, and noticed that their Civil Boundaries ArcMap has linked to it Census population counts for 1990, 2000, 2010 and 2020 for those interested in population changes for counties and localities without having to download and link the Census data. Get it here:
… shows your location as a blue dot in relationship to campsites, trails, state land. This online map requires cell service, won’t work in remote country, so also bring a paper map.
“I could come up with an interesting idea for a blog post, or a could just write a line or two to connect to a government GIS server, and call that a blog post.”
Earlier today I posted a Leaflet Map showing the current Satellite Image from NOAA’s GEOS-East derived product, the GeoColor. It’s nice but to look really attractive you need to increase the brightness and saturation. If you want a 1920×1080 image that looks nice that you can download, you can use this python code. This ARC GIS Image Server is relatively straightforward to work with — supply a bounding box in Web Mercator, size and a few other parameters and it will spit back the image you want.
import requests, io
from PIL import Image, ImageEnhance, ImageDraw, ImageFont
from datetime import datetime
# new york - new england
url = "https://satellitemaps.nesdis.noaa.gov/arcgis/rest/services/MERGED_GeoColor/ImageServer/exportImage/exportImage?bbox=-9011008.3905,4747656.7008,-7345292.6701,6195679.7647&size=1920%2C1080&format=jpgpng&transparent=false&bboxSR=3857&imageSR=3857&f=image"
# full usa
#url = "https://satellitemaps.nesdis.noaa.gov/arcgis/rest/services/MERGED_GeoColor/ImageServer/exportImage/exportImage?bbox=-14401959.1214%2C2974317.6446%2C-7367306.5342%2C6447616.2099&size=1920%2C1080&format=jpgpng&transparent=false&bboxSR=3857&imageSR=3857&f=image&fbclid=IwAR1HvlVi55Bd933eOhM-Cvou2JgPdKAlcp2go4f8B5pqcioAPvcJyfYT0BA"
response = requests.get(url)
image_bytes = io.BytesIO(response.content)
img = PIL.Image.open(image_bytes)
enhancer = ImageEnhance.Color(img)
img = enhancer.enhance(2.5)
enhancer = ImageEnhance.Brightness(img)
img = enhancer.enhance(1.2)
img.save('/tmp/latest_image.jpg')