RSS Git Download  Clone
Raw Blame History
use 5.34.0;
use Feature::Compat::Class;

#================================================================================
class NearestTown; # uses data from GeoNames.org - generated by merge_geonames.pl
#================================================================================

use Data::Printer;
use IO::All -utf8;    # for file handling
use POSIX qw(round);
use List::Util 'uniq';
use Math::Trig qw(great_circle_distance deg2rad);

# ---------------------------
# CONFIG
# ---------------------------
my $cities_file    = 'cities_full.txt';    # output from merge_geonames.pl
my $oceans_file    = 'hydro_features.txt'; # from GeoNames.org
my $min_population = 5000;                 # ignore tiny settlements
my $search_radius  = 1.0;                  # degrees for bounding box

# ---------------------------
# Load cities into memory
# ---------------------------
field $cities = do {
    my @csrc = io($cities_file)->chomp->slurp;
    my @cities;
    for (@csrc) {
        # 51.24389 -1.26154 Overton Hampshire England United Kingdom GB 3392
        my ($lat,$lon,$name,$admin1,$admin2,$cname,$cc,$pop) = split /\t/;
        next unless $pop =~ /^\d+$/ && $pop >= $min_population; # will skip header row
        push @cities, {
            name     => $name,    # city/town name  
            lat      => $lat + 0, # ensure numeric
            lon      => $lon + 0, # ensure numeric
            cc       => $cc,      # abbreviated country code [GB, ES, US etc]
            country  => $cname,
            admin1   => $admin1,  # state/province   
            admin2   => $admin2,  # not useful
            pop      => $pop + 0, # ensure numeric
        };
    }  # p @cities
    \@cities;
};

# -----------------------------
# Load oceans/seas into memory
# -----------------------------
field $oceans = do {
    my @src = io($oceans_file)->chomp->slurp;
    my @hydro;
    for (@src) {
		next unless $_ =~ /^\d+/; # skip header row
        my @f = split /\t/;
        my ($name, $feat_code, $lat, $lon) = @f[1,3,4,5];
        push @hydro, { name => $name, lat => $lat, lon => $lon, code => $feat_code };
    }
    \@hydro;
};

# ---------------------------
# Distance function (Haversine)
# ---------------------------
method haversine_km {
     # convert degrees to radians
    my ($lat1, $lon1, $lat2, $lon2) = map { deg2rad($_) } @_;
    # compute great-circle distance, multiplying by 6371 converts radians → kilometres
    # since 6371 km ≈ Earth’s mean radius.
    return 6371 * great_circle_distance( $lon1, $lat1, $lon2, $lat2 );
}

# ---------------------------
# Nearest town lookup
# ---------------------------
method nearest_town ($lat, $lon) { # p $lat; p $lon;
    # creates a latitude/longitude rectangle centered on your tile ($lat, $lon), 
    # extending $search_radius degrees in each direction. Example: if $lat = 60,
    # $search_radius = 0.5, then the box covers 59.5–60.5°N.   
    my ($min_lat, $max_lat) = ($lat - $search_radius, $lat + $search_radius);
    my ($min_lon, $max_lon) = ($lon - $search_radius, $lon + $search_radius);
    # $best will hold the closest city record found.
    # $min_dist starts at the maximum possible distance within the tile
    my $tile_diagonal_km = 0.5 * sqrt( (111)**2 + (111 * cos(deg2rad($lat)))**2 );         
    my ($best, $min_dist) = (undef, $tile_diagonal_km);

    for my $c (@$cities) {
        # This is a fast bounding-box prefilter — saves time by ignoring cities too 
        # far away (no need to compute distance). 
        next if $c->{lat} < $min_lat || $c->{lat} > $max_lat;
        next if $c->{lon} < $min_lon || $c->{lon} > $max_lon; # p $c->{lon};
        # compute great-circle distance using centre of tile - lat + 0.5, lon + 0.5
        my $d = $self->haversine_km($lat + 0.5, $lon + 0.5,  $c->{lat}, $c->{lon});
           #  printf "Tile diagonal km: %.1f\n", $tile_diagonal_km;
            # print what you think you used
            # printf "Tile centre: %dN %dE\n", $lat + 0.5, $lon + 0.5;

            # print city coordinates and the computed distance:
            # printf "Candidate: %s (%s) at %.2f, %.2f -> dist=%.1f km\n",
            #    $c->{name}, $c->{country}, $c->{lat}, $c->{lon}, $d;
       # If this city is closer than any previous one, store it as the current best match.
        if ( $d < $min_dist ) { $best = $c; $min_dist = $d; } 
    }
    # say "No nearby city found for $lat / $lon" unless $best;
	if ($best) {
		my $formatted_info = $self->format_city_info($best, $min_dist);
		return $formatted_info;
	}
	else {
		my $hydro_name = $self->is_ocean_tile($lat, $lon);
		return "$hydro_name [ocean tile]" if $hydro_name;
	}
    return 0; # can't find anything    
}

method format_city_info ($city, $dist) {
    my $location = join ', ', grep { $_ } uniq( @{$city}{ qw/name admin1 admin2 country/ } ); # not using cc
    return sprintf '%s [%s km from tile center]', $location, round($dist);
}

method is_ocean_tile ($lat, $lon, $threshold_km = 50) { # default 50 km
    for my $o (@$oceans) {
        my $d = $self->haversine_km($lat, $lon, $o->{lat}, $o->{lon});
		if ( $d < $threshold_km ) {
			my $hydro_name = $o->{name};
			# append sea / ocn / bay unless already contains it
			$hydro_name .= ' ' . ucfirst lc $o->{code} 
				unless $hydro_name =~ /$o->{code}\Z/i;
			return $hydro_name;
		}
    }
    return 0; # can't find anything
}

1;