<jupyter_start><jupyter_text>ipyleaflet - NAD27Testing the NAD27 datum with ipyleaflet<jupyter_code>from ipyleaflet import Map, basemaps, basemap_to_tiles, projections, GeoData, DrawControl, LayersControl, FullScreenControl, WMSLayer
import geopandas as gpd
from ipywidgets import Label, HBox, VBox, Layout
from pyproj import Transformer
import matplotlib.pyplot as plt
from sidecar import Sidecar<jupyter_output><empty_output><jupyter_text>The following cell is to avoid putting the data we need in the repo. The idea is that we have multiple datasets and we only downlod what we need for each notebook.<jupyter_code>from catalog import import_data
import_data()
#EPSG 32040 NAD27 Stateplane TX South Central
EPSG_32040 = {
    'name': 'EPSG:32040',
    'custom': True,
    'proj4def': '+proj=lcc +lat_1=28.38333333333333 +lat_2=30.28333333333333 +lat_0=27.83333333333333 +lon_0=-99 +x_0=609601.2192024384 +y_0=0 +datum=NAD27 +units=us-ft +no_defs',
    'origin': [1934715.99110518, 505965.75122003],
    'bounds': [
        [60341.26161686, 30289.16709928],
        [3757196.70267378, 1073405.75104801]
    ],
    'resolutions': [
        8192.0,
        4096.0,
        2048.0,
        1024.0,
        512.0
    ]
}
# So you we what the lat lon would be for regular WGS84 datum
transformer_geographic = Transformer.from_crs(4267, 4326)
# The projected coord in ft from the projection origin
transformer_projected = Transformer.from_crs(4267, 32040)
# This is to see the map on the side panel
sc = Sidecar(title='ipyleaflet map')
coord_nad27_label = Label()
coord_wgs84_label =  Label()
coord_32040_label =  Label()
feature_label = Label()

def handle_interaction(**kwargs):
    if kwargs.get('type') == 'mousemove':
        native_coordinates = kwargs.get('coordinates')
        w84_coordinates = transformer_geographic.transform(native_coordinates[0], native_coordinates[1])
        # Here, leaflet by default will give you the geographic coordinate but you can transform it into your projected coord one just as easily.
        p32040_coordinates = transformer_projected.transform(native_coordinates[0], native_coordinates[1])
        coord_nad27_label.value = f' NAD27 coordinates: {str(native_coordinates)}'
        coord_wgs84_label.value = f' WGS84 coordinates: {str(list(w84_coordinates))}'
        coord_32040_label.value = f' Projected coordinates(ft) 32040: {str(list(p32040_coordinates))}'

m = Map(center=(32, -104),
        zoom=2,
        min_zoom=1,
        layers=(),
        scroll_wheel_zoom=True,
        crs=EPSG_32040)

m.on_interaction(handle_interaction)
m.add_control(LayersControl(position='topright'))
m.add_control(FullScreenControl())

# Looks like there is a small bug in the Sidecar panel
# after using the map full screen and back to normal the bug goes away
with sc:
    display(VBox([m,
                  coord_nad27_label,
                  coord_wgs84_label,
                  coord_32040_label,
                  feature_label],
                  layout=Layout(display='flex',
                                flex_flow='column',
                                border='solid 2px',
                                align_items='stretch',
                                height='100%'))
           )
#COUNTIES_PERMIAN_AREA CRS: NAD27
selDf = gpd.read_file('data/texas-counties/COUNTIES_PERMIAN_AREA.shp')
selDf.head()
selDf.crs
#Convert to NAD27 Stateplane TX South Central
geoDf = selDf.to_crs(epsg=32040)
geoDf.head()
geoDf.crs
fig, axs = plt.subplots(1, 2, figsize=(10, 5))

selDf.plot(ax=axs[0], legend=True)
axs[0].set_title('EPSG:4267')
geoDf.plot(ax=axs[1], legend=True)
axs[1].set_title('EPSG:32040')
plt.show()
def handle_feature_click(feature, **kwargs):
    feature_label.value = str(feature['properties']['CNAME'])

# You don't need to use the reprojected dataframe. Leaflet uses geodetic lat lons and reprojects vector
geoData = GeoData(geo_dataframe=selDf, name='Counties')
geoData.on_click(handle_feature_click)
m.add_layer(geoData)<jupyter_output><empty_output><jupyter_text>Because leaflet was already initialized as 32040 you don't need to reproject your dataframe. It knows how to project it and get you back the proper coordinate using NAD27. There are caveats with different datums, especially if you want to do other operations with "regular" data based on WGS84.The next cell uses a OpenStreet map layer as base layer but their WMS doesn't support 32040 so you'll see the inconsistencies.I couldn't find a good WMS service that supported state plane projections for base layers.<jupyter_code>wms = WMSLayer(url='https://ows.terrestris.de/osm/service',
               name='Unprojected base layer',
               layers='OSM-Overlay-WMS',
               min_zoom=1,
               tile_size=256)
m.add_layer(wms)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0951885023953343, 'fraction_numerical': 0.07081857946261196, 'mean_word_length': 3.7544554455445542, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 13, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 3, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '352618', 'license': 'MIT', 'n_tokens_mistral': 1799, 'n_tokens_neox': 1558, 'n_words': 411, 'path': 'notebooks/nad27test.ipynb', 'repo_name': 'betolink/ipyleaflet-projections', 'size': '4801'}