I have a dataset in the WGS84 geographic coordinate system that I would like to interpolate along a grid using Python.
It appears gdal_grid
is not available through the Python bindings, so I plan to use another 2D interpolation method after rasterizing the data.
To rasterize the data, I plan to re-project it using pyproj, then write it to a raster using the gdal python bindings. The UTM projection coordinate system seems like it would be the best projection, but my data spans the North Atlantic (i.e. multiple UTM zones).
I'm open to suggestion on the interpolation method, but I had planned to use an inverse distance interpolation (as shown here). These are presence absence sightings data from ship cruise tracks, which do not necessarily follow lines of latitude. They span from latitudes in the mid 70's to mid 50's.
What projection would be best to use?
Answer
Sightings constitute a (non-random) sample of some process or population. Accordingly, interpolation (especially) IDW is not a good idea: it solves a different problem altogether.
Consider making a density map. When doing so, it's probably better to favor equal-area projections over conformal projections (because changes of area bias the density, whereas non-conformality does not, even though it might change the qualitative appearance of the map).
There are many equal-area projections of the world. (Unfortunately, popular projections, including Mercator, Transverse Mercator, and Stereographic, are not among them.) Most of them are cylindrical. This is nice in one way--lines of latitude and longitude form an orthogonal grid on the map--but it's bad near the poles, because the scale distortion is all focused in one (vertical) direction where it can become extreme.
These considerations suggest a compromise: just about any projection suitable for a largish area (like a large country), when shifted to the center of the North Atlantic, would likely do a good job, simply because its overall relative distortion should be limited. For instance, the Albers equal-area conic projection is popular for the conterminous US. Changing its central meridian to -30 degrees, its reference parallel to 55 degrees north, and its standard parallels to 45 and 65 degrees, re-points it at the North Atlantic:
The Tissot indicatrices show how little distortion actually occurs within this region: scale errors are less than 1.5% in any direction throughout; that's much more accuracy than needed for interpolation or kernel density estimates.
In contrast--to see what really has been accomplished by this choice compared to others that might suggest themselves--consider the Mercator projection:
The huge discrepancies between the inner and outer circles of the Tissot indicatrices testify to the large changes in scale from south to north in this map. In effect, any interpolation or density calculation based on these map coordinates would use neighborhoods near the top that are approximately half the size of the neighborhoods near the bottom (in reality, not on the map). That's appreciable distortion; there's no sense introducing it into the calculations when it's simple enough to choose a better projection.
No comments:
Post a Comment