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

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.

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 =~/^#/;

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.

• Ok so for me to fully understand. Creating that pass point table I can use any 4 known points on the map? I don't exactly now the borders but some of the circles on the map I know the X|Y as well as Lat|Lon. And since my goal is not to show the map but to use the image and plot additional Lat|Lon on it (the same image). How can I with a map like that re-map additional Lat|Lon to X|Y? Dec 28, 2021 at 18:27
• @fkrauthan OK for clarification, you want to transform something new from lon/lat you know under the discussed parameter set, to the pix.x & pix.y on the image itself and draw some pix based stuff? Dec 28, 2021 at 19:15
• Yes correct. I have this images with a couple of known points. And now I want to draw on that original image a track (for simplicity a straight line is fine) between brand new Lat|Lon coordinates. Dec 28, 2021 at 19:21
• @fkrauthan In which language, JavaScript? Dec 28, 2021 at 19:29
• Yep that would be the best as I am doing it inside the browser. Dec 28, 2021 at 19:45