Dealing with geospatial affine transformations

My goal is to produce a map like the attachment, in which the affine adjustments were done in Postgres/PostGIS, within Julia.


This is a standard cartographic convention for maps of the United States.

The attached map uses the identical shapefile that I cannot successfully transform using various combinations of Shapefile to bring in as typeof polygons, GeoDataFrame to bring in as typeof IGeometry. I’ve tried using GeometryBasics, LibGEOS, ArchGDAL and others. The most common error is that object conversion returns typeof nothing.

For present purposes, I can export by PostGIS derived geometry as GeoJSON to allow users without a PostGIS server to prepare maps, but I would prefer to be able to script it.

It would be helpful to hear the community’s experience.

  1. Is it possible?
  2. Has it been done?
  3. If it has, can you point me to a script?
  4. If there is no public script, is there a workflow known to work

Thanks for any help.

You can do that with GMT. The fig would be made with 3 layers: mainland, Hawaii and Alaska. The trick is to pick the map limits of the Hawaii and Alaska layers such that, when plotted at the same map scale (I imagine you want this), they have the same paper width.

1 Like

An example with a colder Hawaii and a warmer Alaska

using GMT

coast(R="-147/-109/29/51", proj=(name=:lambertConic, center=[-130 40], parallels=[35 45]),
          figsize=15, borders=2, area=500, shore=true)

# Alaska
coast!(R="-167.17/-131.72/53/72.9", figsize=4, borders=2, area=500, shore=true, xshift=2, yshift=1)

# Hawaii
coast!(R="-161/-154/17/24", figsize=2.5, shore=true, xshift=3, yshift=5, show=1)

Thanks. I wanted to stay within Makie so I ended up taking a deep dive into projections and bounding boxes, coming up with this:

using CairoMakie
using GeoDataFrames
using DataFramesMeta
import GeometryOps as GO
using GeoInterface
using Proj

# Define the Albers Equal Area projection for continental US
const AEA_CONUS = "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"

# Define projections for Alaska and Hawaii (centered on their regions)
const AEA_ALASKA = "+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"
const AEA_HAWAII = "+proj=aea +lat_1=8 +lat_2=18 +lat_0=13 +lon_0=-157 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"

# Function to reproject geometry using Proj.jl
function reproject_geom_proj(geom, source_proj4, target_proj4)
    # Create transformation
    trans = Proj.Transformation(source_proj4, target_proj4, always_xy=true)
    
    # Transform the geometry
    return GO.transform(geom) do point
        x, y = trans(point[1], point[2])
        return (x, y)
    end
end

# Improved fit_to_bbox function that works with projected coordinates
function fit_to_bbox_projected(geom, target_min, target_max; rotation_degrees=0)
    # Get geometry bounds
    extent = GeoInterface.extent(geom)
    min_pt = (extent.X[1], extent.Y[1])
    max_pt = (extent.X[2], extent.Y[2])
    
    # Compute scale factors
    scale_x = (target_max[1] - target_min[1]) / (max_pt[1] - min_pt[1])
    scale_y = (target_max[2] - target_min[2]) / (max_pt[2] - min_pt[2])
    scale = min(scale_x, scale_y)  # uniform scaling
    
    # Compute center-based transformation
    center_src = ((min_pt[1] + max_pt[1]) / 2, (min_pt[2] + max_pt[2]) / 2)
    center_tgt = ((target_min[1] + target_max[1]) / 2, (target_min[2] + target_max[2]) / 2)
    
    # Convert rotation to radians
    rotation_rad = rotation_degrees * π / 180
    
    # Apply transformation: translate to origin, scale, rotate, translate to target
    function transform_point(p)
        # Translate to origin and scale
        x, y = (p[1] - center_src[1]) * scale, (p[2] - center_src[2]) * scale
        
        # Apply rotation if specified
        if rotation_degrees != 0
            x_rot = x * cos(rotation_rad) - y * sin(rotation_rad)
            y_rot = x * sin(rotation_rad) + y * cos(rotation_rad)
            x, y = x_rot, y_rot
        end
        
        # Translate to target position
        return (x + center_tgt[1], y + center_tgt[2])
    end
    
    # Apply to all points in geometry
    return GO.transform(transform_point, geom)
end

# Main plotting function
function plot_us_with_insets(states_gdf, alaska_gdf, hawaii_gdf; 
                            zoom_to_insets=false,
                            title="",
                            caption="",
                            legend_elements=nothing,
                            legend_labels=nothing,
                            legend_title="",
                            colors=nothing,
                            show_north_arrow=false,
                            show_scale_bar=false)
    # Separate continental US states
    conus_states = @chain states_gdf begin
        @rsubset(:STUSPS ∉ ["AK", "HI", "PR", "VI", "GU", "MP", "AS"])
    end
    
    # Define source CRS (NAD83 geographic)
    source_crs = "+proj=longlat +datum=NAD83 +no_defs"
    
    # Reproject continental US to Albers Equal Area
    conus_geoms_projected = [reproject_geom_proj(g, source_crs, AEA_CONUS) for g in conus_states.geometry]
    
    # Reproject Alaska
    ak_geom = alaska_gdf.geometry[1]  # assuming single geometry
    # First project Alaska using its own Albers projection
    ak_geom_alaska_aea = reproject_geom_proj(ak_geom, source_crs, AEA_ALASKA)
    # Then reproject to continental US projection for consistent units
    ak_geom_conus_proj = reproject_geom_proj(ak_geom, source_crs, AEA_CONUS)
    
    # Reproject Hawaii
    hi_geom = hawaii_gdf.geometry[1]  # assuming single multipolygon
    
    # Filter to just the main inhabited islands (first 10 polygons)
    # Get the component polygons
    hi_polygons = GI.getgeom(hi_geom)
    
    # Take only the first 10 polygons (the inhabited islands)
    inhabited_polygons = collect(Iterators.take(hi_polygons, 10))
    
    # Create a new MultiPolygon with just the inhabited islands
    hi_geom_filtered = GI.MultiPolygon(inhabited_polygons)
    
    # First project Hawaii using its own Albers projection
    hi_geom_hawaii_aea = reproject_geom_proj(hi_geom_filtered, source_crs, AEA_HAWAII)
    # Then reproject to continental US projection
    hi_geom_conus_proj = reproject_geom_proj(hi_geom_filtered, source_crs, AEA_CONUS)
    
    # Get extent of continental US for positioning
    # Combine all conus geometries to get full extent
    all_conus_points = []
    for geom in conus_geoms_projected
        ext = GeoInterface.extent(geom)
        push!(all_conus_points, (ext.X[1], ext.Y[1]))
        push!(all_conus_points, (ext.X[2], ext.Y[2]))
    end
    
    conus_min_x = minimum(p[1] for p in all_conus_points)
    conus_max_x = maximum(p[1] for p in all_conus_points)
    conus_min_y = minimum(p[2] for p in all_conus_points)
    conus_max_y = maximum(p[2] for p in all_conus_points)
    
    conus_width = conus_max_x - conus_min_x
    conus_height = conus_max_y - conus_min_y
    
    # Calculate a visual reference line for inset placement
    # Instead of using the absolute bottom (Key West), use a line roughly at the 
    # bottom of the southwestern states (approximately 31-32°N)
    # This is roughly 7-8 degrees north of Key West, or about 15-20% up from the bottom
    visual_reference_y = conus_min_y + conus_height * 0.15  # Adjust this percentage as needed
    
    # Debug output
    println("Continental US bounds - X: $conus_min_x to $conus_max_x, Y: $conus_min_y to $conus_max_y")
    println("Visual reference line (SW states bottom): $visual_reference_y")
    println("Difference from true bottom: $(visual_reference_y - conus_min_y) meters")
    
    # Define target bounding boxes (in projected meters)
    # Alaska: even smaller for conventional layout
    alaska_width = conus_width * 0.15  # Smaller than before
    alaska_height = conus_height * 0.15  # Smaller than before
    # Position Alaska further left with tighter gap
    alaska_x_offset = -conus_width * 0.05  # Less left offset to bring closer to Hawaii
    # Use gap from visual reference line (SW states) not from Key West
    alaska_y_gap = 100000  # 100km gap from visual reference line
    alaska_target_min = (conus_min_x + alaska_x_offset, visual_reference_y - alaska_height - alaska_y_gap)
    alaska_target_max = (conus_min_x + alaska_x_offset + alaska_width, visual_reference_y - alaska_y_gap)
    
    # Debug output
    println("Alaska Y gap: $alaska_y_gap")
    println("Alaska height: $alaska_height")
    println("Visual reference Y (SW states): $visual_reference_y")
    println("Alaska top edge: $(visual_reference_y - alaska_y_gap)")
    println("Alaska bottom edge: $(visual_reference_y - alaska_height - alaska_y_gap)")
    println("Distance from visual reference to Alaska top: $(alaska_y_gap)")
    println("Alaska bounds - X: $(alaska_target_min[1]) to $(alaska_target_max[1]), Y: $(alaska_target_min[2]) to $(alaska_target_max[2])")
    
    # Hawaii: make much bigger now that we only have inhabited islands
    hawaii_width = conus_width * 0.45  # Much bigger - was 0.4
    hawaii_height = conus_height * 0.45  # Much bigger
    hawaii_x_center = conus_min_x + conus_width * 0.3  # Position
    hawaii_y_offset = 0  # Keep at same level as Alaska
    hawaii_target_min = (hawaii_x_center - hawaii_width/2, alaska_target_min[2] + hawaii_y_offset)
    hawaii_target_max = (hawaii_x_center + hawaii_width/2, alaska_target_max[2] + hawaii_y_offset)
    
    # Apply transformations to Alaska and Hawaii
    ak_transformed = fit_to_bbox_projected(ak_geom_conus_proj, alaska_target_min, alaska_target_max)
    hi_transformed = fit_to_bbox_projected(hi_geom_conus_proj, hawaii_target_min, hawaii_target_max, rotation_degrees=30)
    
    # Create the plot with layout grid
    fig = Figure(size = (1200, 800), backgroundcolor = :white)
    
    # Main map axis - use subplot grid for better control
    ax = Axis(fig[2, 1:3], aspect = DataAspect(), backgroundcolor = :lightblue)
    
    # Plot continental US
    if isnothing(colors)
        # Default colors if none provided
        for geom in conus_geoms_projected
            poly!(ax, geom, color = :white, strokecolor = :black, strokewidth = 0.5)
        end
    else
        # Use provided colors for continental US
        conus_colors = colors[1:length(conus_geoms_projected)]
        for (i, geom) in enumerate(conus_geoms_projected)
            poly!(ax, geom, color = conus_colors[i], strokecolor = :black, strokewidth = 0.5)
        end
    end
    
    # Plot Alaska with appropriate color
    ak_color = isnothing(colors) ? :white : colors[length(conus_geoms_projected) + 1]
    poly!(ax, ak_transformed, color = ak_color, strokecolor = :black, strokewidth = 0.5)
    
    # Plot Hawaii with appropriate color
    hi_color = isnothing(colors) ? :white : colors[length(conus_geoms_projected) + 2]
    poly!(ax, hi_transformed, color = hi_color, strokecolor = :black, strokewidth = 0.5)
    
    # Add bounding boxes for insets (optional)
    alaska_bbox_points = [
        Point2(alaska_target_min...),
        Point2(alaska_target_max[1], alaska_target_min[2]),
        Point2(alaska_target_max...),
        Point2(alaska_target_min[1], alaska_target_max[2]),
        Point2(alaska_target_min...)
    ]
    
    hawaii_bbox_points = [
        Point2(hawaii_target_min...),
        Point2(hawaii_target_max[1], hawaii_target_min[2]),
        Point2(hawaii_target_max...),
        Point2(hawaii_target_min[1], hawaii_target_max[2]),
        Point2(hawaii_target_min...)
    ]
    
    lines!(ax, alaska_bbox_points, color = :transparent, linewidth = 1, linestyle = :dash)
    lines!(ax, hawaii_bbox_points, color = :transparent, linewidth = 1, linestyle = :dash)
    
    # Add labels
    text!(ax, "                           Alaska and Hawaii not to scale", position = (alaska_target_min[1] + alaska_width/2, alaska_target_max[2] + 50000),
          align = (:center, :bottom), fontsize = 14, font = "Arial Bold")
    text!(ax, "", position = (hawaii_x_center, hawaii_target_max[2] + 50000),
          align = (:center, :bottom), fontsize = 14, font = "Arial Bold")
    
    # Set axis limits based on zoom preference
    if zoom_to_insets
        # Zoom to show just the bottom portion with insets
        plot_min_y = min(alaska_target_min[2], hawaii_target_min[2]) - 50000
        plot_max_y = conus_min_y + conus_height * 0.3  # Show only bottom 30% of continental US
        println("Zoomed view: focusing on insets and southern states")
    else
        # Show full map
        plot_min_y = min(alaska_target_min[2], hawaii_target_min[2]) - 50000
        plot_max_y = conus_max_y + 50000
    end
    
    plot_min_x = min(alaska_target_min[1], conus_min_x) - 50000
    plot_max_x = max(hawaii_target_max[1], conus_max_x) + 50000
    
    xlims!(ax, plot_min_x, plot_max_x)
    ylims!(ax, plot_min_y, plot_max_y)
    
    # Debug output
    println("Plot Y limits: $plot_min_y to $plot_max_y (range: $(plot_max_y - plot_min_y))")
    println("Plot X limits: $plot_min_x to $plot_max_x (range: $(plot_max_x - plot_min_x))")
    println("Continental US height: $conus_height")
    println("Y-axis spans $(round((plot_max_y - plot_min_y)/1000)) km")
    println("Expected Y-span should be roughly $(round((conus_height + alaska_height + alaska_y_gap + 100000)/1000)) km")
    
    # Remove axis decorations
    hidedecorations!(ax)
    hidespines!(ax)
    
    # Add title if provided
    if !isempty(title)
        Label(fig[1, :], title, fontsize = 24, font = "Arial Bold")
    end
    
    # Add legend if elements provided
    if !isnothing(legend_elements) && !isnothing(legend_labels)
        Legend(fig[3, 3], legend_elements, legend_labels, legend_title;
               halign = :right, valign = :top, fontsize = 12)
    end
    
    # Add caption if provided
    if !isempty(caption)
        Label(fig[4, :], caption, fontsize = 10, halign = :right)
    end
    
    # Add north arrow if requested
    if show_north_arrow
        text!(fig.scene, "⇧\nN", 
              position=(0.85, 0.6),
              space=:relative,
              align=(:left, :top), 
              fontsize=36, color=:black)
    end
    
    # Add scale bar if requested
    if show_scale_bar
        # Calculate scale bar for 100 miles
        map_width_meters = conus_width  # width in meters
        miles_to_meters = 160934.4  # 100 miles in meters
        scale_bar_length = (miles_to_meters / map_width_meters) * 0.8  # as fraction of plot width
        
        scale_bar_x = 0.85
        scale_bar_y = 0.45
        
        lines!(fig.scene, 
               [scale_bar_x - scale_bar_length/2, scale_bar_x + scale_bar_length/2], 
               [scale_bar_y, scale_bar_y], 
               space=:relative, color=:black, linewidth=3)
        
        text!(fig.scene, "100 mi", 
              position=(scale_bar_x, scale_bar_y - 0.02),
              space=:relative,
              align=(:center, :top), 
              fontsize=14, color=:black)
    end
    
    return fig
end

# Example usage:
# Assuming you have loaded your shapefiles into GeoDataFrames
states = GeoDataFrames.read("data/2024_shp/cb_2024_us_state_500k.shp")
alaska = @rsubset(states, :STUSPS == "AK")
hawaii = @rsubset(states, :STUSPS == "HI")

fig = plot_us_with_insets(states, alaska, hawaii, title="Map of the United States with Alaska and Hawaii insets")
#save("us_map_with_insets.png", fig)
![us_map_with_insets|690x459](upload://jbIXBtDbkfZEPAB2UH4g98kqucA.png)

Since this is for general chloropleth mapping, I decided not to include the graticule

GeoMakie supports a lot of transformations:

Not sure if it will fit your exact use case though :wink:

Certainly true. I think I went through most of them. :sweat_smile: I ended up with this because of several challenges, including getting the right bounding box for Hawaii (it’s state shapefile goes all the way to Midway). Anyway, it was fun.

1 Like

@technocrat yet another alternative in native Julia is to load the shapefiles with GeoIO.jl and call the Proj and Affine transforms from GeoStats.jl to perform the transformations you need. It should be much more intuitive than handling all these C PROJ strings manually:

using GeoStats
using GeoIO

import GLMakie as Mke

shp1 = GeoIO.load("file1.shp")
shp2 = GeoIO.load("file2.shp")

res1 = shp1 |> Proj(Albers) |> Affine(A1, b1)
res2 = shp2 |> Proj(Albers) |> Affine(A2, b2)

viz(res1.geometry)
viz!(res2.geometry)

I will leave the details up to you.

1 Like

Thank you. I had a bit of a hike to construct the projections for Alaska and Hawaii, but that eventuated in

const AlaskaAlbers = Albers{50, 55, 65, NAD83}
projector_ak = Proj(AlaskaAlbers)
alaska |> projector_ak 

Since it’s relatively to find CRS strings in the form

"+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"

once I understood the use of the constructor, I was able to get there.

This is, as you note, far better than what I had attempted. Thanks, again.

1 Like

You may need an extra shift for the lon_0 parameter:

crs = CoordRefSystems.shift(Albers{50,55,65}, lonₒ=-154)
1 Like

@technocrat notice that this coordinate reference system corresponds to EPSG code 3338:

I’ve added the code to CoordRefSystems.jl here:

It should be available soon as a patch release:

After you update your environment, you will also be able to use the code in the Proj transform:

shp |> Proj(EPSG{3338})

Anytime you find a code that is not mapped yet, please feel free to submit a similar PR.

2 Likes