The dataset is documented in "MANAIR Manual of Standards and Procedures for Aviation Forecasts: appendix C" -- Appendix C: Graphic Area Forecast (GFA)—Supplement to Chapter 4. A coordinate system is not mentioned here but it seem to be a polar stereo-graphic projection.
Lets give it a try:
1. Find ground gontrol points (GCP's) and create a pass point table
The border polygon can be used to set ground control points (GCP's) in the graphic (at least four) (from GADM for example),
Create a pass point table (QGIS has a geo-referencing tool for example.)
(Corner points are used ..purple marked circles)
Initial longitude/latitude based pass point table:
# BC Table GCP
# NAME PIX.X PIX.Y WORLD.LON WORLD.LAT
GCP1 99 90 -139 60
GCP2 368 137 -120.05 60
GCP3 485 453 -114.06 49
GCP4 312 459 -122.75 49
2. Calculate corresponding pass points at the spherical coordinate system:
Apply the proj.4 definition string of to get the GCP corresponding to the projection.
The proj.4 string is:
proj +proj=stere +lat_0=90 +lat_ts=60
-139.00 60.00<CR>
-2097489.33 2412885.46
...
The final table looks like that (..see Perl script at the end):
PIX: 99 90 LON.LAT: -139.00 60.00 SPHERE: -2097489.33 2412885.46
PIX: 368 137 LON.LAT: -120.05 60.00 SPHERE: -2767377.74 1600967.90
PIX: 485 453 LON.LAT: -114.06 49.00 SPHERE: -4070468.40 1817399.23
PIX: 312 459 LON.LAT: -122.75 49.00 SPHERE: -3749152.78 2411536.11
3. Geo-referencing with gdal_translate:
Try to geo-reference your picture by these points using the tool gdal_translate
# SET SRS COMMAND:
/usr/bin/gdal_translate \
-gcp 099 090 -2097489.33 2412885.46 \
-gcp 368 137 -2767377.74 1600967.89 \
-gcp 485 453 -4070468.40 1817399.23 \
-gcp 312 459 -3749152.78 2411536.10 \
-a_srs '+proj=stere +lat_0=90 +lat_ts=60' \
-expand rgb \
-of GTiff \
data/Latest-gfacn31_cldwx_000.png \
data/Latest-gfacn31_cldwx_000_trans.tif
with:
-gcp -- ground control point: img.col img.row world.x world.y,-a_srs -- define the coordinate system (proj.4 string),-expand -- expand from a palette define color image to RGB 24 bit image,-of --output format (refer to gdal formats for further info's),data/Latest-gfacn31_cldwx_000.png -- incoming dataset anddata/Latest-gfacn31_cldwx_000_trans.png -- outgoing dataset.
4. Resampling to have a good image quality with gdalwarp:
Resample the image to an defined resolution and coordinate system here Canadian Polar Stereographic (EPSG:5937) with gdalwarp
# RESAMPLE COMMAND:
/usr/bin/gdalwarp \
-overwrite \
-t_srs '+init=epsg:5937' \
-tr 2000 2000 \
-r lanczos \
data/Latest-gfacn31_cldwx_000.tif_trans.tif \
data/Latest-gfacn31_cldwx_000.tif-warp.tif
with:
-t_srs -- tagret coordinate system (proj.4 string again),-tr 2000 2000 -- target resolution [m], here 2x2 km,-r lancoz -- Interpolation method here sinc(x) based,data/Latest-gfacn31_cldwx_000.tif_trans.tif -- incoming dataset anddata/Latest-gfacn31_cldwx_000.tif-warp.tif.
4. GIS Application qgis:
Now you can use your "Aviation Weather Map" in an Geo-Information-System like QGIS, overlay spatial data and create augmented maps.
5. Perl Sketch
I've written a small Perl script to get a more compact workflow and used the tools directly by a shell call (..same like bash $(...)).
There are similar tools for JavaScript proj4js or [npm gdal(https://www.npmjs.com/package/gdal)] or in QGIS more functional available.
The script and the related datasets can be found under https://github.com/bigopensky/gis-se-snippet.
#!/usr/bin/env perl
use strict;
use warnings;
use Geo::Proj4;
use Data::Dumper;
# Find the tools in the operation system
my $toolTrans = `which gdal_translate`;
my $toolWarp = `which gdalwarp`;
chomp($toolTrans);
chomp($toolWarp);
# Incoming file
my $inFile = 'data/Latest-gfacn31_cldwx_000.png';
# Temporary file
my $outTemp = $inFile;
# Strip extension
$outTemp =~ s/\.png$//;
# Tanslate step to set the srs
my $outTrans = "${outTemp}_trans.tif";
# Warp step to oversample the stuff
my $outWarp = "${outTemp}_warp.tif";
# Source coordinate system
my $srcCrs='+proj=stere +lat_0=90 +lat_ts=60';
# Final coordinate system
my $dstCrs='+init=epsg:5937';
# Destination resolution
my $gsd =2000; # 2 km
# Set the psspoint input as text
my $GCP_LL='
# ID PY WLON WLAT
GPC1 099 090 -139.00 60.00
GPC2 368 137 -120.05 60.00
GPC3 485 453 -114.06 49.00
GPC4 312 459 -122.75 49.00';
# Create the projection
my $proj = Geo::Proj4->new($srcCrs)
or die "parameter error: ".Geo::Proj4->error. "\n";
# Read the coordinate system lines
my @coords = split(/\n/, $GCP_LL);
my @gcps;
for my $coord (@coords) {
# Trim data string
$coord =~ s/^\s+//;
$coord =~ s/\s+$//;
# Next line if empty
next if $coord eq '';
# Next line if comment
next if $coord =~/^#/;
# Read coordinates
my ($id, $px, $py, $wlon, $wlat) = split(/\s+/, $coord);
# Transform lon lat to sphere x and y
my ($wx, $wy) = $proj->forward($wlat, $wlon);
printf "PIX: %3d %3d LON.LAT: %7.2f %6.2f SHP: %2.2f %2.2f\n",
$px, $py, $wlon, $wlat, $wx, $wy;
push (@gcps,"-gcp $px $py $wx $wy");
}
# Translate GCP int parameter strings
my $gcpstr = join(" ", @gcps);
my ($cmd, $res);
# Set the reference system and ground control points
$cmd ="$toolTrans $gcpstr -a_srs '$srcCrs' -expand rgb -of GTiff $inFile $outTrans";
print "SET SRS COMMAND: $cmd\n";
$res = `$cmd`;
# Resample the image
$cmd = "$toolWarp -overwrite -t_srs '$dstCrs' -tr $gsd $gsd -r lanczos $outTrans $outWarp" ;
print "RESAMPLE COMMAND: $cmd\n";
$res = `$cmd`;
# EOF
EDIT: Comment Browser based JavaScript Approach of the Original
If the original image should be used as a "Map Canvas" for a browser based application, the product with source CRS and affine transform data produced by gdal_translate (data/Latest-gfacn31_cldwx_000.tif_trans.tif in the scripts above) can be applied in conjuction with libraries/applications like Leaflet to draw the Lon/Lat based vector content.
Examples: https://github.com/stuartmatthews/leaflet-geotiff
