Lambert Conformal Conic to Geographic Transformation Formulae Land Information New Zealand
Publication No. 47: Lambert conformal conic projection with two standard parallels d/l as QB275U35no47.pdf
Publication No. 49: Lambert projection tables with conversion tables (supplement to 47, above, d/l as QB275U35no491918.pdf)
Publication No. 53: General Theory of the Lambert Conformal Conic Projection Oscar Adams, 1918
Publication No. 52, Lambert projection tables for the United States (1 item available at NOAA Miami Library.)
Publication No. 57: General Theory of Polyconic Projections (d/l as: general_theory_polygonic_57.pdf)
penn.edu The Online Books Page
.../filedvsflown/vuelayers_v12/src/views
SINCE we can't install onto rserver without approval, we'll fire up our own venv:
$ cd openmaptiles/
$ cd localpython/
$ source bin/activate
(localpython) $ pip3.7 install pycrs
(localpython) $ pip3.7 install pyproj==2.6.1
https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1
>>> from pyproj import Transformer
>>> transformer = Transformer.from_crs("epsg:4326", lcc, always_xy=True)
>>> transformer.transform( -97.0376947222222, 32.8972330555556)
(136861.77805443306, -121274.43604351794)
SOME help from here:
custom map projection demo #447
https://github.com/ghettovoice/vuelayers/issues/447
/home/wturner/ciws/meekma/meow/etl/sectionals
the clean_proj.sh script sed's out the projection coords from the html text
views/reproject.py is script to transform coords
SEARCHing epsg.org:
this one seems close, but I don't want it:
Name: NAD27 / US National Atlas Equal Area
Code: 9311
CRS Type: Projected
Name: NAD83
Code: 4269
CRS Type: Geographic 2D
Scope: Geodesy.
Extent: North America - NAD83Open
SEARCHing for USA gets 1900 results
SEARCHing for USA + LCC gets 12 results for KS, OR, VA
SEARCHing for USA + Lambert gets 13 results for TX, VA, OR
there is an esri one defined in pyproj, but not in EPSG, so other systems (using proj) don't recognize it
$ python3.7
>>> from pyproj import CRS
>>> cc = CRS('esri:102004')
>>> cc
<Projected CRS: ESRI:102004>
Name: USA_Contiguous_Lambert_Conformal_Conic
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: USA - CONUS - onshore
- bounds: (-124.79, 24.41, -66.91, 49.38)
Coordinate Operation:
- name: USA_Contiguous_Lambert_Conformal_Conic
- method: Lambert Conic Conformal (2SP)
Datum: North American Datum 1983
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich
those limits are:
HOWEVER, epsg.io DOES list this:
http://epsg.io/102004
ESRI:102004
USA Contiguous Lambert Conformal Conic
Description: +proj=lcc +lat_1=33 +lat_2=45 +lat_0=39 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs
the WKT as html shows:
PROJCS["USA_Contiguous_Lambert_Conformal_Conic",
GEOGCS["GCS_North_American_1983",
DATUM["North_American_Datum_1983",
SPHEROID["GRS_1980",6378137,298.257222101]],
PRIMEM["Greenwich",0],
UNIT["Degree",0.017453292519943295]],
PROJECTION["Lambert_Conformal_Conic_2SP"],
PARAMETER["False_Easting",0],
PARAMETER["False_Northing",0],
PARAMETER["Central_Meridian",-96],
PARAMETER["Standard_Parallel_1",33],
PARAMETER["Standard_Parallel_2",45],
PARAMETER["Latitude_Of_Origin",39],
UNIT["Meter",1],
AUTHORITY["EPSG","102004"]]
which seems ok, as does PROJ.4:
+proj=lcc
+lat_1=33
+lat_2=45
+lat_0=39
+lon_0=-96
+x_0=0
+y_0=0
+datum=NAD83
+units=m
+no_defs
https://desktop.arcgis.com/en/arcmap/10.3/guide-books/map-projections/lambert-conformal-conic.htm
says:
Limitations
Best for regions predominantly east - west in extent and located in the middle north or south latitudes. Total latitude range should not exceed 35 d.
Uses and applications SPCS for all zones that are more east - west in extent.
USGS 71/2 - minute quad sheets to match the State Plane Coordinate System.
Used for many new USGS maps created after 1957. It replaced the Polyconic projection.
Continental United States: standard parallels, 33 d and 45 d N.
donwloaded:
Map Projections - A Working Manual USGS
p 104 - 110 describes LCC
Invalid Projection EPSG:102008
https://github.com/geopandas/geopandas/issues/943
ESRI:102004
USA Contiguous Lambert Conformal Conic
http://epsg.io/102004
OpenLayers transform from docs
https://openlayers.org/en/latest/apidoc/module-ol_proj.html#.transform
Using coordinates in meters on a static image map
https://github.com/ghettovoice/vuelayers/issues/399
which is:
var oldCode = OSM_layer.getSource().getProjection().getCode(); // Getting projection of target layer
newCoord = ol.proj.transformExtent([extent.j[0], extent.j[1], extent.j[2], extent.j[3]], oldCode, 'EPSG:4326'); // perform projection transform with given extent.
stackexchange example of transform:
https://gis.stackexchange.com/questions/238243/transform-coordinates-epsg32628-to-epsg3857-openlayers-3
Transforming coordinates of feature in openlayers https://stackoverflow.com/questions/36134974/transforming-coordinates-of-feature-in-openlayers
INSERTING my own definition:
ESRI:102004
USA Contiguous Lambert Conformal Conic
srid: 102004
auth_name: WT's LCC (ESRI)
auth_srid: 102004
srtext:
PROJCS["USA_Contiguous_Lambert_Conformal_Conic",
GEOGCS["GCS_North_American_1983",
DATUM["North_American_Datum_1983",
SPHEROID["GRS_1980", 6378137, 298.257222101]],
PRIMEM["Greenwich", 0],
UNIT["Degree",0.017453292519943295]],
PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["False_Easting",0],
PARAMETER["False_Northing",0],
PARAMETER["Central_Meridian",-96],
PARAMETER["Standard_Parallel_1",33],
PARAMETER["Standard_Parallel_2",45],
PARAMETER["Latitude_Of_Origin",39],
UNIT["Meter",1],
AUTHORITY["EPSG","102004"]
]
proj4text: +proj=lcc +lat_1=33 +lat_2=45 +lat_0=39 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs
and to actually insert into db=naturalearth table=spatial_ref_sys on opal4
and into db=meekma table=spatial_ref_sys on rserver:
INSERT INTO spatial_ref_sys VALUES(
102004,
'WTs LCC (ESRI)',
102004,
'PROJCS["USA_Contiguous_Lambert_Conformal_Conic",GEOGCS["GCS_North_American_1983",DATUM["North_American_Datum_1983",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",-96],PARAMETER["Standard_Parallel_1",33],PARAMETER["Standard_Parallel_2",45],PARAMETER["Latitude_Of_Origin",39],UNIT["Meter",1],AUTHORITY["EPSG","102004"]]',
'+proj=lcc +lat_1=33 +lat_2=45 +lat_0=39 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs');
INSERT 0 1
ok, our own LCC is inserted into postgis on my own opal4 and rserver
formulas:
https://gis.stackexchange.com/questions/53940/calculating-bounding-box-from-known-centre-coordinate-and-zoom
https://community.jmp.com/t5/Discussions/Is-it-possible-to-generate-a-bounding-box-around-a-set-of/td-p/51031
lame: https://www.mapranksearch.com/
which says:
Given you're looking for a simple formula, this is probably the simplest way to do it, assuming that the Earth is a sphere with a circumference of 40075 km.
Length in meters of 1 d of latitude = always 111.32 km
Length in meters of 1 d of longitude = 40075 km * cos( latitude ) / 360
good: https://projections.mgis.psu.edu/
Name Lambert Conic Conformal (2SP)
EPSG Code 9802
GeoTIFF Code CT_LambertConfConic_2SP (8)
CT_LambertConfConic (8)
OGC WKT Name Lambert_Conformal_Conic_2SP
Supported By EPSG, GeoTIFF, PROJ.4, OGC WKT
python pgm to convert: https://gis.stackexchange.com/questions/144061/convert-coordinates-to-differnt-spatial-reference-system-in-python
and this one:
it seems that the LCC bounding box is measured from the central meridian, i.e. +lon_0=-96 and the "central parallel", which must be half-way between lat_1 and lat_2, which is 39.
THUS, a very tiny bounding box would be:
lower left lonlat: POINT (-96.0013 38.999 0)
lower left lcc : POINT (-112.001118326269 -110.408614161706 0)
upper right lonlat: POINT (-95.9987 39.001 0)
upper right lcc : POINT (111.997959401995 110.410232447696 0)
so, it seems that .001 in the Y direction is about equal to .0013 in the X direction
NOTE: these LCC meter coords are very close to a square (centered around coords of 0,0):
coords from def in pyproj for esri:102004:
then, using:
https://www.openstreetmap.org/export#map=4/31.24/-86.66
on rserver:
cd openmaptiles/localpython/
source bin/activate
cd ../lambertconformal/
./transf.py
and using that a LOT, we finally arrived at this as a bounding box:
enter ll: -2367543 -2555495
lcc : POINT (-2367543 -2555495 0)
lonlat: POINT (-116.369826129393 14.1495659747164 0)
enter ur: 2310000 2122048
lcc : POINT (2310000 2122048 0)
lonlat: POINT (-61.0494711865487 53.9467981319802 0)
width: 4677543.0
height: 4677543.0
squareness: 0.0
put this into pg_tileserv/config/ilc_config.toml:
Xmin = -2367543
Ymin = -2555495
Xmax = 2310000
Ymax = 2122048
put this into OL_nascnrs.vue on rserver:
proj_my_own_lcc.setExtent([ -2367543, -2555495, 2310000, 2122048 ]);
1) AFTER inserting 102004 into rserver's spatial_ref_sys table, we can do this:
SELECT ST_AsText(ST_Transform(wkb_geometry,102004)) from atac_sectors limit 1;
POLYGON((-829272.879176356 -209464.414991489,-799043.743155599 -212525.478097802,...
2) and we edit draw_high.py and it seems fine.
3) and to convert each of centers, ultra, low:
ogr2ogr -f "PostgreSQL" PG:"dbname='meekma'" \
-s_srs EPSG:4326 \
-t_srs "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=39 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs" \
/home/data/shapefiles/ARTCC_boundaries/ffs_only_usdom.shp \
-nln "centers_lcc"
worked fine, but so what? Why not do the ST_Transform on the read???
ALSO, the column def is wrong:
meekma=# \d centers_lcc
Table "public.centers_lcc"
Column | Type | Collation | Nullable |
--------------+--------------------------+-----------+----------+
ogc_fid | integer | | not null |
area | numeric(12,3) | | |
perimeter | numeric(12,3) | | |
name | character varying(20) | | |
nas | character varying(12) | | |
id | numeric(3,0) | | |
minalt | numeric(16,4) | | |
maxalt | numeric(16,4) | | |
wkb_geometry | geometry(Polygon,900914) | | |
Indexes:
"centers_lcc_pkey" PRIMARY KEY, btree (ogc_fid)
"centers_lcc_wkb_geometry_geom_idx" gist (wkb_geometry)
Q: convert all, and have a .gjsn for 4326 AND a .gjsn for 102004, or put them all into a postgis table, select, st_translate as needed, and ST_AsGeoJSON via sql???
FOR NOW, we'll create xxx_lcc.gjson files and see how well/bad that works out...
EDITS for lcc:
filedvsflown/analysis/geojson/all_shp_to_geojson.sh
filedvsflown/etl/utils/atac/draw_high.py
Will attempt to make LCC bounaries larger:
approx like this:
> pip3 install ogr
> pip3 install wheel
> pip3 install osr
Successfully built osr
Installing collected packages: Pyglet, coverage, deprecation, dataclasses, filetype, eyeD3, osr
Attempting uninstall: dataclasses
Found existing installation: dataclasses 0.7
Uninstalling dataclasses-0.7:
Successfully uninstalled dataclasses-0.7
betterproto 2.0.0b3 requires dataclasses<0.8,>=0.7, but you'll have dataclasses 0.8 which is incompatible.
Successfully installed Pyglet-1.1.4 coverage-5.5 dataclasses-0.8 deprecation-2.1.0 eyeD3-0.9.6 filetype-1.0.8 osr-0.0.1
/home/wendell/openmaptiles/vmpython/bin/python -m pip install --upgrade pip
NOPE, was too hard to install (osr, osgeo, gdal, etc)
transf.py on rserverthis is looking good:
enter ll: -3400000 -3400000
lcc : POINT (-3400000 -3400000 0)
lonlat: POINT (-122.717813367545 5.44183678260296 0)
enter ur: 3100000 3100000
lcc : POINT (3100000 3100000 0)
lonlat: POINT (-43.271122202429 58.2627068165145 0)
width: 6500000.0
height: 6500000.0
squareness: 0.0
Xmin = -3400000
Ymin = -3400000
Xmax = 3100000
Ymax = 3100000
proj_my_own_lcc.setExtent([ -3400000, -3400000, 3100000, 3100000 ]);
Q: why didn't you figure out that it was square, just like Paul said???
anyway, here is the new one:
