1

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.

1 Answer 1

2

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.)

enter image description here

(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.

enter image description here

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

Comparence GIS vs. WEB Canvas

10
  • 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?
    – fkrauthan
    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?
    – huckfinn
    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.
    – fkrauthan
    Dec 28, 2021 at 19:21
  • @fkrauthan In which language, JavaScript?
    – huckfinn
    Dec 28, 2021 at 19:29
  • Yep that would be the best as I am doing it inside the browser.
    – fkrauthan
    Dec 28, 2021 at 19:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.