Map lat/lon onto polar stereographic projection map based on two known points on map image

Question

I've searched already everywhere but couldn't find anything. I have a "polar stereographic projection" map that is true at 60 degrees latitude north as an image (e.g. https://flightplanning.navcanada.ca/Latest/gfa/anglais/produits/uprair/gfa/gfacn31/Latest-gfacn31_cldwx_000.png). I know on that image the X/Y pixel coordinates as well as the corresponding Lat/Lon coordinates for a couple of points (but not the actual edges).

Now I am trying to figure out how I can map a couple of Lat/Lon coordinates onto that map. So I need to convert Lat/Lon to X/Y pixel.

Currently trying todo that in JavaScript inside the browser to visualize a path on a chart like that.

Answer 1

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 and
• data/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 and
• data/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