https://github.com/spacetelescope/spherical_geometry
https://github.com/anutkk/sphericalgeometry
good definition and explanation: Lambert Conformal Conic
lecture notes: Choosing a projection
ArcGIS: Choose the right projection
wikipedia: Azimuthal equidistant projection
Map Projections Used by the U.S. Geological Survey, 1982
Algorithms for geodesics Charles F. F. Karney; Journal of Geodesy volume 87, pages 43–55(2013)
Geodesics on an Ellipsoid web resources for Karney paper.
WILL USE: Lambert Azimuthal Equal Area
someone said: To map middle latitudes, use a conic projection
Use geographiclib for accurate geodesic distance and area calcuations
https://code.env.duke.edu/projects/mget/ticket/549
https://mrjean1.github.io/PyGeodesy/
https://github.com/mrJean1/PyGeodesy
PyGeodesy A pure Python implementation of geodesy tools for various ellipsoidal and spherical earth models using precision trigonometric, vector-based, elliptic and approximate methods for geodetic (lat-/longitude) and geocentric (ECEF cartesian) coordinates.
duke.edu: Use geographiclib for accurate geodesic distance and area calcuations includes examples of distance(?)
s.o.:
Looking for a pythonic way to calculate the length of a WKT linestring
says:
The geopy module provides the Vincenty formula, which provides accurate ellipsoid distances. Couple this with the wkt loading in Shapely, and you have reasonably simple code:
geopy/geopy is mostly a geocoding service, but as a "by-product" it has a geodesice distance function included:
>>> from geopy.distance import great_circle
>>> newport_ri = (41.49008, -71.312796)
>>> cleveland_oh = (41.499498, -81.695391)
>>> print(great_circle(newport_ri, cleveland_oh).miles)
536.997990696
so, it seems to only work on single segements.
PyGeodesy A pure Python implementation of geodesy tools for various ellipsoidal and spherical earth models using precision trigonometric, vector-based, elliptic and approximate methods for geodetic (lat-/longitude) and geocentric (ECEF cartesian) coordinates.
Package pygeodesy detailed function description
Find the intersection between two geographical data points
geographiclib This implements Algorithms for Geodesics (Karney, 2013) for solving the direct and inverse problems for an ellipsoid of revolution.
GeographicLib initial C++ lib
Transforming Shapely Polygon and MultiPolygon objects
so: How to calculate the area of a polygon on the earth's surface using python? An interesting discussion on spherical vs elipsoid, and an interesting conversion, yet uses manual conversion(???)
list of stuff: Awesome Geospatial Awesome Long list of geospatial analysis tools. Geospatial analysis, or just spatial analysis, is an approach to applying statistical analysis and other analytic techniques to data which has a geographical or spatial aspect.
not sure, probably lame: python-geospatial
not sure, probably lame: spatialpandas Pandas and Dask extensions for vectorized spatial and geometric operations.
list of libraries: Open Source Spatial Analysis Tools for Python: A Quick Guide
list of companies: Awesome Geospatial Companies
another (lame) list of libraries: Essential geospatial Python libraries
pyogrio - Vectorized spatial vector file format I/O using GDAL/OGR
use ogr2ogr to convert shapefile:
ogr2ogr -t_srs EPSG:4326 out.shp in.shp
ISSUES:
Lambert Conic Conformal (2SP)
seems to indicate that 9802 is a generic 2SP LCC code,
and apparently requires
the Central Meridian and the two Standard Parallels to work properly.
that means that the proj.4 string will be:
+proj=lcc +lat_1=38.6666666666667
+lat_2=33.3333333333333
+lat_0=37.5
+lon_0=-107.5
+x_0=0
+y_0=0
BASEGEOCRS:
CONVERSION:
$ gdalinfo Denver SEC 104.tif
Driver: GTiff/GeoTIFF
Files: Denver SEC 104.tif
Size is 16635, 12391
Coordinate System is:
PROJCRS["Lambert Conformal Conic",
BASEGEOGCRS["NAD83",
DATUM["North American Datum 1983",
ELLIPSOID["GRS 1980",6378137,298.257222101004,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4269]],
CONVERSION["Lambert Conic Conformal (2SP)",
METHOD["Lambert Conic Conformal (2SP)",
ID["EPSG",9802]],
PARAMETER["Latitude of false origin",37.5,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8821]],
PARAMETER["Longitude of false origin",-107.5,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8822]],
PARAMETER["Latitude of 1st standard parallel",38.6666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8823]],
PARAMETER["Latitude of 2nd standard parallel",33.3333333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8824]],
PARAMETER["Easting at false origin",0,
LENGTHUNIT["metre",1],
ID["EPSG",8826]],
PARAMETER["Northing at false origin",0,
LENGTHUNIT["metre",1],
ID["EPSG",8827]]],
CS[Cartesian,2],
AXIS["easting",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["northing",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]
Data axis to CRS axis mapping: 1,2
Origin = (-380521.035404817666858,291529.854634510295000)
Pixel Size = (42.334666405531003,-42.334857259625188)
Metadata:
AREA_OR_POINT=Area
TIFFTAG_DATETIME=2020:12:04 10:29:04
TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
TIFFTAG_SOFTWARE=Adobe Photoshop CC 2019 (Windows)
TIFFTAG_XRESOLUTION=300
TIFFTAG_YRESOLUTION=300
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( -380521.035, 291529.855) (111d57'15.54"W, 40d 2'50.37"N)
Lower Left ( -380521.035, -233041.362) (111d41'25.51"W, 35d19'26.36"N)
Upper Right ( 323716.140, 291529.855) (103d42'35.65"W, 40d 4' 8.14"N)
Lower Right ( 323716.140, -233041.362) (103d56' 4.30"W, 35d20'39.77"N)
Center ( -28402.448, 29244.246) (107d49'21.13"W, 37d45'47.67"N)
Band 1 Block=16635x1 Type=Byte, ColorInterp=Palette
Color Table (RGB with 256 entries)
0: 255,255,255,255
1: 255,255,0,255
2: 255,0,255,255
...
253: 104,104,104,255
254: 24,24,24,255
255: 8,8,8,255
sqlite3 /home/data/local/share/proj/proj.db "SELECT code FROM projected_crs WHERE auth_name = 'EPSG';"