Geo Python

python
geopandas
Author

Geo

Published

November 2, 2023

OK, first off, something OhhelloAna said inspired me to remember something…

I will never forget solving a difficult problem that had stumped me for days, in just a few moments after a relaxing three-day weekend 🍸😺 You have inspired me to write about it on my blog & I will do it as soon as I think of a good picture🍸🐈

So here is a picture of a restful scene, to remind everyone to take time off regularly to rest, to play, to re-randomize, and of course to hydrate. At the beach.

restful scene of life-giving water

I don’t remember the thing I solved back then, but it involved putting a lot of data onto a CD-ROM. So … it was tough, innovative, and quite a long time ago. Computers don’t even have DVD drives in theme anymore.

Onward to minimum viable python land!

Hey Bing, can you …

… write a python program to convert a geographic shapefile to EPSG-4326?

Bing says “Sure, here’s a Python script that uses the pyproj and geopandas libraries to accomplish this. Please make sure to install these libraries using pip before running the script.”

I fire up Pycharm, community edition, stumble around getting the code into a project and getting the new python version into the environment, and voila. PyCharm says there are errors and I click on each of them to let it help out by installing geopandas, etc. Later, I improve the script to take command line arguments, and Bing is helpful in telling how to configure PyCharm CE to help debug those, and then that script works too. It’s a manual for everything!

import geopandas as gpd

from pyproj import CRS

# Load the shapefile

gdf = gpd.read_file('ShapeIn/B.shp')

# Check the current coordinate system

print(f"Original Coordinate System: {gdf.crs}")

# Convert the geometry to EPSG-4326

–> THIS IS THE LINE <–

gdf = gdf.to_crs(CRS("EPSG:4326"))

# Check the new coordinate system

print(f"New Coordinate System: {gdf.crs}")

# Write the modified shapefile

gdf.to_file('ShapeOut/B.shp')

So, here is a basically a python one-liner … the pgm comes with ‘insert your path here’, which is what I do because i have a file called B.shp in the “ShapeIn” folder … and it totally works. I mean, you could quibble with the line gdf = gdf.to_crs(CRS("EPSG:4326"))

but we’re beggars, right? Thanks, ChatGTP4

I did this because I had written an R program to do the same thing, but for some odd reason R Studio is doing bad things on Windows10 … a white screen of death that they believe they have fixed … so I thought I would explore the python options, if I didnt’ always want to work in MacOS. (VERY hypothetical. VERY.)

Hey Anaconda, you got a real one-liner?

Install anaconda

Anaconda fails to download to Windows, so I go to the mac.

only for you, etc. etc. Running package scripts, etc etc.

Then, with conda resident, I install odg

install odg (gdal.org)

conda install -c conda-forge gdal

https://casoilresource.lawr.ucdavis.edu/software/gdal-and-ogr-geodata-conversion-and-re-projection-tools/

I read the Davis website that shows how to operate ogr2ogr, and there are many fun things, like with -append, you can append to an existing shapefile, and that you can just specify the transform with -t_srs EPSG:4326, so that’s great. In fact that is why I have downloaded all these gigabytes, for that one line.

So I first try out ogrinfo, and it spits out a humongous list of lat-long coordinates and much meta info.

>ogrinfo -al ShapeIn/B.shp

Then, grokking the weird order, where out is before in, I also note you can make a KML file from an shapefile, or vice versa. And it will try to guess the input and output types from the filename extension, so yay!

>ogr2ogr -f”KML” outputfilename, inputfilename

>ogr2ogr -f”ESRI shapefile” outputfilename, inputfilename

Anyway, now I am thinking that after installing Anaconda, and using Conda to install gdal from Conda Forge, I have got a one-liner

ogr2ogr -t_srs EPSG:4326 ShapeOut/B_modified.shp ShapeIn/B.shp

and try that, and I get

“Warning 1: One or several characters couldn’t be converted correctly from UTF-8 to ISO-8859-1. This warning will not be emitted anymore.”

BUT

my output shapefile has been created, and it loads successfully.

The operation of re-projection has been performed, which means that ogr2ogr’s helpful message about how there has to be an existing coordinate system for it to be transformed has been fulfilled.

(A source SRS must be available for reprojection to occur.)

And so i have gotten down to a one-liner and if that isn’t minimum viable python, I don’t know what is.

— all photos Copyright © 2022-2024 George D Girton all rights reserved