1package authoring23import (4 "fmt"5 "math"6 "net/url"7 "strconv"8 "strings"910 "gopkg.in/yaml.v3"11)1213const (14 usTopoService = "https://basemap.nationalmap.gov/arcgis/rest/services/USGSTopo/MapServer/export"15 tigerRoads = "https://tigerweb.geo.census.gov/arcgis/rest/services/TIGERweb/Transportation/MapServer"16 tigerHydro = "https://tigerweb.geo.census.gov/arcgis/rest/services/TIGERweb/Hydro/MapServer"17 tigerCounties = "https://tigerweb.geo.census.gov/arcgis/rest/services/TIGERweb/State_County/MapServer/1"18 webMercatorR = 20037508.34278924419 webMercatorLat = 85.0511287820 maxExportSize = int64(4096)21)2223// AreaProfile combines an arbitrary geographic selection with the built-in24// official U.S. source pack. Bounds are WGS84 west, south, east, north25// coordinates; they are not restricted to the pack's coverage.26type AreaProfile struct {27 ID string28 Title string29 Bounds [4]float6430 DetailZoom int31 IncludeTopo bool32 IncludeRoads bool33 IncludeHydro bool34 IncludeCounties bool35 Presentation AreaPresentation36}3738// AreaPresentation is the creator-owned layer, label, ordering, and color39// choice that is compiled into the Atlas presentation and legend.40type AreaPresentation struct {41 RoadLabel, HydroLabel, CountyLabel string42 RoadColor, HydroColor, HydroFill, CountyColor string43 RoadVisible, HydroVisible, CountyVisible bool44 RoadOrder, HydroOrder, CountyOrder int6445}4647// AreaSummary is the cost preview shown before a manifest is installed.48type AreaSummary struct {49 OriginX int6450 OriginY int6451 SourceZoom int52 PixelSize int6453 RasterTiles int6454 CaptureRequests int55 EstimatedBytes int6456}5758// Deprecated compatibility aliases. New code should name the geographic59// selection separately from the U.S.-coverage source pack.60type USAreaProfile = AreaProfile61type USAreaPresentation = AreaPresentation62type USAreaSummary = AreaSummary6364// DefaultAreaProfile selects the built-in official U.S.-coverage source pack.65func DefaultAreaProfile() AreaProfile {66 return AreaProfile{67 ID: "sample-region", Title: "Sample Region",68 Bounds: [4]float64{-105.1, 39.6, -104.9, 39.8},69 DetailZoom: 12, IncludeTopo: true, IncludeRoads: true,70 IncludeHydro: true, IncludeCounties: true,71 Presentation: AreaPresentation{72 RoadLabel: "Roads", HydroLabel: "Water", CountyLabel: "Counties",73 RoadColor: "#b45f3c", HydroColor: "#4b8db8", HydroFill: "#b9d9ea", CountyColor: "#766b5b",74 RoadVisible: true, HydroVisible: true, CountyVisible: true,75 RoadOrder: 40, HydroOrder: 20, CountyOrder: 10,76 },77 }78}7980// DefaultUSAreaProfile is kept for callers that name the bundled source pack.81// Deprecated: use DefaultAreaProfile.82func DefaultUSAreaProfile() USAreaProfile { return DefaultAreaProfile() }8384// NewAreaProject turns one arbitrary Web Mercator region selection into a85// complete portable project using the selected source families. Captures and86// artifacts remain outside the manifest.87func NewAreaProject(profile AreaProfile) (Project, AreaSummary, error) {88 window, err := usTileWindow(profile.Bounds, profile.DetailZoom)89 if err != nil {90 return Project{}, AreaSummary{}, err91 }92 if !profile.IncludeTopo && !profile.IncludeRoads && !profile.IncludeHydro && !profile.IncludeCounties {93 return Project{}, AreaSummary{}, fmt.Errorf("select at least one source")94 }95 if err := validateUSPresentation(profile); err != nil {96 return Project{}, AreaSummary{}, err97 }98 project := usProjectBase(window, profile.DetailZoom, profile.ID, profile.Title)99 if profile.IncludeTopo {100 project.Rasters = append(project.Rasters, usTopoRaster(window, profile.DetailZoom))101 }102 if profile.IncludeRoads {103 addRoads(&project, window.transform, profile.Presentation)104 }105 if profile.IncludeHydro {106 addHydro(&project, window.transform, profile.Presentation)107 }108 if profile.IncludeCounties {109 addCounties(&project, window.transform, profile.Presentation)110 }111 if err := project.Validate(); err != nil {112 return Project{}, AreaSummary{}, fmt.Errorf("validate area project: %w", err)113 }114 return project, summarizeUSProject(project, window), nil115}116117// NewUSAreaProject is kept for callers that name the bundled source pack.118// Deprecated: use NewAreaProject.119func NewUSAreaProject(profile USAreaProfile) (Project, USAreaSummary, error) {120 return NewAreaProject(profile)121}122123func validateUSPresentation(profile AreaProfile) error {124 choices := []struct {125 enabled bool126 label string127 colors []string128 order int64129 }{130 {profile.IncludeRoads, profile.Presentation.RoadLabel, []string{profile.Presentation.RoadColor}, profile.Presentation.RoadOrder},131 {profile.IncludeHydro, profile.Presentation.HydroLabel, []string{profile.Presentation.HydroColor, profile.Presentation.HydroFill}, profile.Presentation.HydroOrder},132 {profile.IncludeCounties, profile.Presentation.CountyLabel, []string{profile.Presentation.CountyColor}, profile.Presentation.CountyOrder},133 }134 for _, choice := range choices {135 if !choice.enabled {136 continue137 }138 if strings.TrimSpace(choice.label) == "" || choice.order < 0 {139 return fmt.Errorf("enabled presentation layers require labels and non-negative order")140 }141 for _, color := range choice.colors {142 if !validHexColor(color) {143 return fmt.Errorf("presentation color %q must be #RRGGBB", color)144 }145 }146 }147 return nil148}149150func validHexColor(value string) bool {151 if len(value) != 7 || value[0] != '#' {152 return false153 }154 for _, character := range value[1:] {155 if !((character >= '0' && character <= '9') || (character >= 'a' && character <= 'f') || (character >= 'A' && character <= 'F')) {156 return false157 }158 }159 return true160}161162// MarshalProjectYAML emits the single portable manifest after validating it.163func MarshalProjectYAML(project Project) ([]byte, error) {164 if err := project.Validate(); err != nil {165 return nil, err166 }167 data, err := yaml.Marshal(project)168 if err != nil {169 return nil, fmt.Errorf("encode project manifest: %w", err)170 }171 return data, nil172}173174type usWindow struct {175 originX, originY int64176 side, pixels int64177 maxZoom int64178 bbox [4]float64179 transform CoordinateTransform180}181182func usTileWindow(bounds [4]float64, zoom int) (usWindow, error) {183 for _, coordinate := range bounds {184 if math.IsNaN(coordinate) || math.IsInf(coordinate, 0) {185 return usWindow{}, fmt.Errorf("area coordinates must be finite")186 }187 }188 if bounds[0] >= bounds[2] || bounds[1] >= bounds[3] {189 return usWindow{}, fmt.Errorf("area requires west < east and south < north")190 }191 if bounds[0] < -180 || bounds[2] > 180 || bounds[1] < -webMercatorLat || bounds[3] > webMercatorLat {192 return usWindow{}, fmt.Errorf("area must stay within the Web Mercator world")193 }194 if zoom < 6 || zoom > 16 {195 return usWindow{}, fmt.Errorf("detail zoom must be between 6 and 16")196 }197 n := int64(1) << zoom198 minX, maxY := lonTile(bounds[0], n), latTile(bounds[1], n)199 maxX, minY := lonTileUpper(bounds[2], n), latTileUpper(bounds[3], n)200 originX, originY, side := alignedSquare(minX, minY, maxX, maxY)201 pixels := side * 256202 if pixels > maxExportSize {203 return usWindow{}, fmt.Errorf("selection needs a %d-pixel source image; maximum is %d (lower detail or narrow the area)", pixels, maxExportSize)204 }205 bbox := mercatorTileBounds(originX, originY, side, n)206 scale := float64(256*n) / (2 * webMercatorR)207 transform := CoordinateTransform{Kind: "affine", Matrix: [6]float64{208 scale, 0, webMercatorR*scale - float64(256*originX),209 0, -scale, webMercatorR*scale - float64(256*originY),210 }}211 return usWindow{originX: originX, originY: originY, side: side, pixels: pixels, maxZoom: int64(math.Log2(float64(side))), bbox: bbox, transform: transform}, nil212}213214func lonTile(longitude float64, count int64) int64 {215 return int64(math.Floor((longitude + 180) / 360 * float64(count)))216}217218func lonTileUpper(longitude float64, count int64) int64 {219 return int64(math.Ceil((longitude+180)/360*float64(count))) - 1220}221222func latTile(latitude float64, count int64) int64 {223 value := (1 - math.Asinh(math.Tan(latitude*math.Pi/180))/math.Pi) / 2 * float64(count)224 return int64(math.Floor(value))225}226227func latTileUpper(latitude float64, count int64) int64 {228 value := (1 - math.Asinh(math.Tan(latitude*math.Pi/180))/math.Pi) / 2 * float64(count)229 return int64(math.Floor(value))230}231232func alignedSquare(minX, minY, maxX, maxY int64) (int64, int64, int64) {233 for side := int64(1); ; side *= 2 {234 originX, originY := minX-(minX%side), minY-(minY%side)235 if originX+side-1 >= maxX && originY+side-1 >= maxY {236 return originX, originY, side237 }238 }239}240241func mercatorTileBounds(originX, originY, side, count int64) [4]float64 {242 unit := 2 * webMercatorR / float64(count)243 return [4]float64{244 -webMercatorR + float64(originX)*unit,245 webMercatorR - float64(originY+side)*unit,246 -webMercatorR + float64(originX+side)*unit,247 webMercatorR - float64(originY)*unit,248 }249}250251func usProjectBase(window usWindow, sourceZoom int, id, title string) Project {252 if id == "" {253 id = "sample-region"254 }255 if title == "" {256 title = "Sample Region"257 }258 return Project{259 Schema: ProjectSchema, SchemaNamespace: "atlas.local/" + id + "/v1",260 ID: id, Title: title,261 Target: Target{World: id, Title: title, CoordinateSpace: CoordinateSpace{262 ID: id + "-grid", Kind: "projected", Unit: "pixel", Definition: "atlas:tile-plane",263 Extent: [4]float64{0, 0, float64(window.pixels), float64(window.pixels)},264 SourceZoom: int64(sourceZoom), OriginX: window.originX, OriginY: window.originY,265 TileSize: 256, Size: window.pixels,266 }},267 Presentation: Presentation{ID: "default", Title: title},268 Release: Release{Revision: 1},269 Budgets: BuildBudgets{RequestBytes: 512 << 20, TotalBytes: 4 << 30, Requests: 10_000, RasterTiles: 10_000, RasterPixels: 1_000_000_000},270 }271}272273func usTopoRaster(window usWindow, sourceZoom int) Raster {274 query := url.Values{275 "bbox": {floatList(window.bbox)}, "bboxSR": {"3857"}, "imageSR": {"3857"},276 "size": {strconv.FormatInt(window.pixels, 10) + "," + strconv.FormatInt(window.pixels, 10)},277 "format": {"png32"}, "transparent": {"false"}, "f": {"image"},278 }279 return Raster{280 ID: "usgs-topo", Name: "USGS Topo", Adapter: "raster-file", Locator: usTopoService + "?" + query.Encode(),281 MediaType: "image/png", License: "U.S. government work; public domain", Attribution: "Map services and data available from U.S. Geological Survey, National Geospatial Program.",282 EstimateBytes: window.pixels * window.pixels, TileSize: 256, SourceZoom: int64(sourceZoom), MaxZoom: window.maxZoom,283 FullZoom: window.maxZoom, Interpolate: true, Background: "#e9e4d7",284 Bounds: &RasterRect{Width: float64(window.pixels), Height: float64(window.pixels)},285 }286}287288func floatList(values [4]float64) string {289 return strconv.FormatFloat(values[0], 'f', 6, 64) + "," + strconv.FormatFloat(values[1], 'f', 6, 64) + "," +290 strconv.FormatFloat(values[2], 'f', 6, 64) + "," + strconv.FormatFloat(values[3], 'f', 6, 64)291}292293func addRoads(project *Project, transform CoordinateTransform, presentation AreaPresentation) {294 project.FeatureSets = append(project.FeatureSets, FeatureSetContract{295 ID: "roads", Title: "Roads", SemanticType: "transportation", Geometry: "path",296 Properties: []PropertyContract{{ID: "name", Name: "Name", Type: "string", Optional: true}, {ID: "class", Name: "Class", Type: "string", Optional: true}, {ID: "route-type", Name: "Route type", Type: "string", Optional: true}},297 })298 for _, layer := range []struct{ id, number string }{{"primary-roads", "2"}, {"secondary-roads", "6"}, {"local-roads", "8"}} {299 source := tigerSource(layer.id, tigerRoads+"/"+layer.number, "roads", "path", transform)300 source.Mapping.Identity = "properties.OID"301 source.Mapping.FeatureTitle = "coalesce(properties.NAME, properties.OID)"302 source.Mapping.Properties = append(source.Mapping.Properties, PropertyMap{Field: "route-type", Source: "properties.RTTYP"})303 project.Sources = append(project.Sources, source)304 }305 project.Presentation.Styles = append(project.Presentation.Styles, Style{ID: "roads", Stroke: presentation.RoadColor})306 project.Presentation.Layers = append(project.Presentation.Layers, Layer{ID: "roads", FeatureSet: "roads", Style: "roads", Label: presentation.RoadLabel, Group: "Transportation", LabelPolicy: "quiet", Visible: presentation.RoadVisible, Order: presentation.RoadOrder, MaxZoom: 32})307}308309func addHydro(project *Project, transform CoordinateTransform, presentation AreaPresentation) {310 project.FeatureSets = append(project.FeatureSets,311 FeatureSetContract{ID: "waterways", Title: "Waterways", SemanticType: "hydrography", Geometry: "path", Properties: commonTigerProperties()},312 FeatureSetContract{ID: "water-bodies", Title: "Water bodies", SemanticType: "hydrography", Geometry: "area", Properties: append(commonTigerProperties(), PropertyContract{ID: "water-area", Name: "Water area", Type: "int64", Optional: true})},313 )314 waterways := tigerSource("waterways", tigerHydro+"/0", "waterways", "path", transform)315 waterBodies := tigerSource("water-bodies", tigerHydro+"/1", "water-bodies", "area", transform)316 waterBodies.Mapping.Properties = append(waterBodies.Mapping.Properties, PropertyMap{Field: "water-area", Source: "properties.AREAWATER"})317 project.Sources = append(project.Sources, waterways, waterBodies)318 project.Presentation.Styles = append(project.Presentation.Styles, Style{ID: "waterways", Stroke: presentation.HydroColor}, Style{ID: "water-bodies", Stroke: presentation.HydroColor, Fill: presentation.HydroFill})319 project.Presentation.Layers = append(project.Presentation.Layers,320 Layer{ID: "water-bodies", FeatureSet: "water-bodies", Style: "water-bodies", Label: presentation.HydroLabel + " bodies", Group: presentation.HydroLabel, LabelPolicy: "quiet", Visible: presentation.HydroVisible, Order: presentation.HydroOrder, MaxZoom: 32},321 Layer{ID: "waterways", FeatureSet: "waterways", Style: "waterways", Label: presentation.HydroLabel + "ways", Group: presentation.HydroLabel, LabelPolicy: "quiet", Visible: presentation.HydroVisible, Order: presentation.HydroOrder + 1, MaxZoom: 32},322 )323}324325func addCounties(project *Project, transform CoordinateTransform, presentation AreaPresentation) {326 project.FeatureSets = append(project.FeatureSets, FeatureSetContract{327 ID: "jurisdictions", Title: "Counties", SemanticType: "jurisdiction", Geometry: "area",328 Properties: []PropertyContract{{ID: "name", Name: "Name", Type: "string"}, {ID: "land-area", Name: "Land area", Type: "int64", Optional: true}, {ID: "water-area", Name: "Water area", Type: "int64", Optional: true}},329 })330 source := tigerSource("counties", tigerCounties, "jurisdictions", "area", transform)331 source.Mapping.Identity, source.Mapping.FeatureTitle = "properties.GEOID", "properties.NAME"332 source.Mapping.Properties = []PropertyMap{{Field: "name", Source: "properties.NAME"}, {Field: "land-area", Source: "properties.AREALAND"}, {Field: "water-area", Source: "properties.AREAWATER"}}333 project.Sources = append(project.Sources, source)334 project.Presentation.Styles = append(project.Presentation.Styles, Style{ID: "counties", Stroke: presentation.CountyColor})335 project.Presentation.Layers = append(project.Presentation.Layers, Layer{ID: "counties", FeatureSet: "jurisdictions", Style: "counties", Label: presentation.CountyLabel, Group: "Context", LabelPolicy: "quiet", Visible: presentation.CountyVisible, Order: presentation.CountyOrder, MaxZoom: 32})336}337338func commonTigerProperties() []PropertyContract {339 return []PropertyContract{{ID: "name", Name: "Name", Type: "string", Optional: true}, {ID: "class", Name: "Class", Type: "string", Optional: true}}340}341342func tigerSource(id, locator, set, family string, transform CoordinateTransform) Source {343 return Source{344 ID: id, Adapter: "arcgis-feature-service", Locator: locator, MediaType: "application/geo+json",345 License: "U.S. government work; public domain", Attribution: "Source: U.S. Census Bureau TIGERweb. This product does not imply Census Bureau endorsement.",346 EstimateBytes: 16 << 20, Query: Query{Where: "1=1", Limit: 1000, ObjectID: "OBJECTID", SpatialReference: 3857},347 Mapping: Mapping{348 FeatureSet: set, Identity: "properties.OBJECTID", FeatureTitle: "coalesce(properties.NAME, properties.OBJECTID)",349 Geometry: GeometryMap{Family: family, SourceSpace: "EPSG:3857", Transform: transform},350 Properties: []PropertyMap{{Field: "name", Source: "properties.NAME"}, {Field: "class", Source: "properties.MTFCC"}},351 },352 }353}354355func summarizeUSProject(project Project, window usWindow) AreaSummary {356 tiles := (int64(1)<<(2*(window.maxZoom+1)) - 1) / 3357 estimated := int64(0)358 for _, source := range project.Sources {359 estimated += source.EstimateBytes360 }361 for _, raster := range project.Rasters {362 estimated += raster.EstimateBytes363 }364 return AreaSummary{365 OriginX: window.originX, OriginY: window.originY, SourceZoom: int(project.Target.CoordinateSpace.SourceZoom),366 PixelSize: window.pixels, RasterTiles: tiles, CaptureRequests: len(project.Sources) + len(project.Rasters), EstimatedBytes: estimated,367 }368}