2

Environment Canada helpfully provides recent captures of all national weather radar stations at http://dd.meteo.gc.ca/radar/PRECIPET/GIF/ at ten minute intervals. These are the same GIF images as those posted on their weather website (e.g. http://weather.gc.ca/radar/index_e.html?id=xsm) and inserted below.

enter image description here

This question usefully describes how to download current weather radar data in KML formats from Environment Canada. However, I need to automate the processing of a deep archive of GIF images captured at ten minute intervals.

I understand, too, that raw radar data in spatial formats is available from Environment Canada at cost recovery pricing that will have much higher thematic and spatial resolution than is possible from these GIF rasters. However, such precision is unimportant for this application.

The challenge then is to define the coordinate system of these GIFs in a known projection. The location of the radar at the centre of the map ("+") is known and published (51.0246, -113.3981 for the radar station in the above image). I have also determined through inquiries to Environment Canada that these radar images are in an azimuthal equidistant projection.

How would one form a PROJ4 string, or define a projection in ArcGIS, that would describe these GIF images in their untransformed state, and so take advantage of this free resource.

Community
  • 1
digitalmaps
  • 501
  • 3
  • 14

2 Answers2

2

ArcGIS supports spheres or ellipsoids (spheroids) with the azimuthal equidistant projection.

Here's a sample prj (wkt) string for "World_Azimuthal_Equidistant" based on WGS 1984.

PROJCS["World_Azimuthal_Equidistant",
GEOGCS["GCS_WGS_1984",
DATUM["D_WGS_1984",
SPHEROID["WGS_1984",6378137.0,298.257223563]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]],
PROJECTION["Azimuthal_Equidistant"],
PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",0.0],
PARAMETER["Latitude_Of_Origin",0.0],
UNIT["Meter",1.0]]

To specify a sphere, you would set the second numerical parameter in the SPHEROID keyword to 0 (zero). As AndreJ added in comments, you would change the central meridian and latitude of origin parameter values to match a particular radar site.

If the gifs all represent the same area in azimuthal equidistant, the corner coordinates (which you'll need to georeference the gifs) will all be the same, no matter what the center point's coordinates are.

mkennedy
  • 18,916
  • 2
  • 38
  • 60
  • The parameters for Central_Meridian and Latitude_Of_Origin should be adjusted to the point given in the question. And for the use with PROJ.4/GDAL/QGIS, a sphere should be used instead of an ellipsoid. – AndreJ Jul 28 '15 at 18:33
  • @AndreJ Good comments! I tried to incorporate them into the answer. – mkennedy Jul 28 '15 at 18:52
  • This is very useful for defining the projection @mkennedy. What would be an efficient way to determine the coordinates of the corners of the GIF, i.e. to georeference the GIF. All I have is the scale bar on the right and,the longlat of the centre point. Features (i.e lakes) with known locations are poorly depicted. Presumably the dimensionality in km of each raster cell could be derived from the scale bar, and then the corners georeferenced using the Euclidean distance in km converted from the raster distance. Is there a more efficient way with ArcGIS? – digitalmaps Jul 28 '15 at 19:38
  • Nothing I know about. You might be able to get NR Canada to give you the general extents. – mkennedy Jul 28 '15 at 19:42
  • 1
    A critical and very useful piece of information is evident upon zooming into the 40 km scale bar on the right. Each pixel in this GIF is square representing 1 km in dimension. It is therefore trivial to determine the corners of the map for georeferencing. The centre is at (240,240) in the GIF coordinate system, so each corner is therefore 339.4 km distant from the known longlat coordinates at the centre of the map. – digitalmaps Jul 28 '15 at 19:59
  • Assuming this map is true in its alignment, then the coordinates of the corners can be found by adding or subtracting 240 km of latitude or longitude given the centre point. – digitalmaps Jul 28 '15 at 20:06
  • @Digital Unfortunately that is incorrect--your prescription assumes lat-lon is a Euclidean coordinate system, which it is not. However, one can use the map to compute the distances of the corners from the center (use the Pythagorean theorem, because all distances from the center are true) as well as their bearings (because the projection is conformal at the center) and then use geodetic forward calculations to find their coordinates. – whuber Jul 29 '15 at 01:27
  • @whuber thank you for this useful caution. Why then would my graphical intuition be incorrect: examining the map above there have been 6 concentric circles spaced at 40 km added by the cartographer. The sixth from the centre is tangent with the map borders, implying that the borders of the map are 240 km from its centre. Understanding that a degree of latitude and longitude are not equidistant, why is not sufficient to separately find the latitude that is 240 km north of the centre point, and the longitude that is 240 km west of the centre point in order to find the top left corner? – digitalmaps Jul 29 '15 at 03:12
  • You can reproject the corner coordinates in the aeqd projection back to lat/long, if it helps. – AndreJ Jul 29 '15 at 04:18
  • @digital Because latitude and longitude do not enjoy the same Pythagorean relationship as x and y in a Euclidean plane. In particular, the top and bottom sides of the map are not lines of latitude. One way to appreciate what might fail is to consider the polar aspect. If the Prime Meridian were straight to the right in the map, then the longitude 240 km west of the center would be 180 degrees, but the longitudes of the left corners would be near +135 and -135 degrees. The discrepancy for oblique aspects, such as the one shown here, is not quite as great but is still substantial. – whuber Jul 29 '15 at 11:33
1

You have to specify the projection on a sphere, because proj seems to only support the spherical formulas of this projection:

+proj=aeqd  +R=6371000 +lat_0=51.0246 +lon_0=-113.3981

For some example plots, see Manipulating Azimuthal Equidistant Projections in QGIS

Community
  • 1
AndreJ
  • 76,698
  • 5
  • 86
  • 162