Geographical Transformations

[kmuehlbauer]

This is a follow up issue on a discussion on the gitter-channel (this comment and ff.) started by @WilfredJanssen about geographical transformations.

While there exists an example about image transformations there is no example which covers geographical transformations.

The only reference to geographical transformations is within vispy/glsl/transforms, but the defined functions are never used in VisPy codebase.

@davidh-ssec did some work building PROJ.4 compatible transform objects in the SIFT project, see some action here.

@WilfredJanssen aims to create a GUI which will be used as some sort of meteorological work station: visualising weather observations (station obs and soundings), satellite and radar data, marine observations, numerical model data etc. etc.

@kmuehlbauer uses VisPY to plot 2D radar data in polar coordinates (azimuth, range) using a polar image transform. But there is also need to plot 3D Volume data in polar stereographic projection or rotated pole projections (like COSMO uses). He uses GDAL/OSR for projection/transformation on the CPU.

This issue should act as the main reference point in the discussion of implementation of geographical transformations in VisPy.

[rougier]

Those transforms come from glumpy where they are actually used in some examples. See for example the choropleth example.

[kmuehlbauer]

@rougier Yes, I knew there was this chloropleth example also in VisPy, but it is not with the geographical transformation. Thanks for the pointer.

[WilfredJanssen]

I already created a function for plotting data on a map which is plotted in a regular lat/lon projection (by far the easiest one).

Data can be plotted by just giving the data and the upper left and lower right corner coordinates. The function itself will place and scale the generated image into the right spot.

In my case, the code is:

plotData(windSpeed, ucLon=-180, ucLat=90, lrLon=180, lrLat=-90, cmap='summer')

The next step is to create several common transformations (like web mercator, lambert etc.). Transforming vector based visuals (e.g. LineVisual) is quite simple since you can calculate the positions with given coordinates. Transforming images (including gridded data) is more complex and requires OpenGL knowledge (the first thing I have to work on).

[djhoese]

In my SIFT project I have a Visual based on LineVisual for reading shapefiles and storing the vertexes as lon/lat coordinates. I have a subclass of the ImageVisual that breaks a large satellite image in to smaller tiles and and different resolutions. It uses PROJ.4 projection strings to define the projection as well as upper-left corner coordinates and size of each pixel in X and Y directions. This is similar to how geotiffs describe their data.

I think @rougier had originally pointed me at some of the glumpy projection stuff, but I didn't know it already existed in the vispy repository. I created my own based on python strings because I needed to be able to customize the projection parameters (which I parse from a PROJ.4 string) in the GLSL. If we can make it easy to do this with loading from glsl files then I like that. I'll have to look more at how glumpy uses them.

So first step forward is to decide how to implement these. Without looking at glumpy I'm guessing there was some reason its implementation wasn't copied over to vispy or maybe it was lower priority so it never made it over. I'm a little biased because this is how I implemented my stuff, but is there any reason that these shouldn't be implemented as a transform object similar to MatrixTransform or STTransform? The forward mapping in these geographic transforms would be lon/lat to projection space, and inverse is...inverse. I had made some arguably hacky choices in my implementation to handle transform errors where I return infinity which has always been ignored in my experience but I'm not sure that is standard for OpenGL.

Lastly, PROJ.4 versus other way of describing things? I'm ok having individual parameters as parameters to subclasses and I could even write a special PROJ4Transform that maps a PROJ.4 string to the individual classes. My implementation in SIFT tries to reduce some of the GPU computations by precomputing some of the calculations in python in the CPU before passing it over to the GLSL.

Brain dump...complete.

[rougier]

@WilfredJanssen For the (very nice by the way) image projection above, do you use the forward or inverse transform? (i.e. do you "raytrace" the image or do you really project it ?)

[WilfredJanssen]

@rougier I project (I guess?) the image using MatrixTransform (for rotating the image 180 degrees) and STTansform for scaling and translating.

The code:

        # Rotate the image 180 degrees
        mapTransformation = MatrixTransform()
        mapTransformation.rotate(180, (1, 0, 0))

        # Image scaling
        xScale = (lowerrightlon - upperleftlon) / xDimension
        yScale = (upperleftlat - lowerrightlat) / yDimension
        # Image translation
        xTranslate = upperleftlon
        yTranslate = (upperleftlat * -1.) + (upperleftlat - lowerrightlat)
        transformation = (STTransform(scale=(xScale, yScale),  translate=(xTranslate, yTranslate, 0)) * 
                          mapTransformation)
        # Set visual transformation
        visual.transform = transformation

[rougier]

I'm asking because after looking at your image, it might make sense to use the inverse transform. If you iterate over all the pixels of the screen image, you might ask what is the corresponding color. It might help speed up the process. I remember having problem when using the forward projection with the apparition of "holes" at place where the distortion was big. But maybe it's a stupid idea (never tested in fact).

[luigi1973]

I tried to implement a transform class for transverse mercator projection. I have simply derived the Vispy BaseTransform class, introducing the projection code into the _gslmap and _gslimap strings. It's a work in progress but it works quite well. Following the code, may be it can be useful for someone.

class TransverseMercator(BaseTransform):
    """ data (degrees,degrees)

    step: default (1,1)
    origin : 
    """
    def __init__(self, origin=(0.0,0.0), step=(1.0, 1.0)):
        super(TransverseMercator, self).__init__()
        self._coords = np.array(origin)/180.0*3.1415
        self._step   = np.array(step)/180.*3.1415
        self._update_shaders()

    glsl_map = """
    vec4 transform_forward(vec4 pos)
    {
    float lambda = (pos.x * $dlon ) + $lon0;
    float phi    = (pos.y * $dlat ) + $lat0;
    float x = 0.5*0.75*log((1.0+sin(lambda)*cos(phi)) / (1.0 - sin(lambda)*cos(phi)));
    float y = 0.75*1.00*atan(tan(phi), cos(lambda));
    return vec4(x,y, pos.z, 1);
    }
    """
    glsl_imap = """
    vec4 transform_inverse(vec4 pos)
    {
    float x = pos.x;
    float y = pos.y;
    float a =1.00;
    float k0 = 0.75;
    float lambda = atan( (0.5 * (exp(x/(k0*a))-exp(-x/(k0*a)))),cos(y/(k0*a)));
    float phi    = asin(sin(y/(k0*a))/ (0.5 * (exp(x/(k0*a))+exp(-x/(k0*a)))));
    return vec4((lambda -$lon0) /$dlon, (phi -$lat0) /$dlat, pos.z+10,1);
    }
    """
    Linear = False
    Orthogonal = False
    NonScaling = False
    Isometric = False

    def _update_shaders(self):
        self._shader_map['lon0'] = self._coords[0]
        self._shader_map['lat0'] = self._coords[1]
        self._shader_map['dlon'] = self._step[0]
        self._shader_map['dlat'] = self._step[1]
        self._shader_imap['lon0'] = self._coords[0]
        self._shader_imap['lat0'] = self._coords[1]
        self._shader_imap['dlon'] = self._step[0]
        self._shader_imap['dlat'] = self._step[1]

[djhoese]

@luigi1973 Thanks. Technically I think you also need to implement the map and imap python functions. Either way, thanks for the code.

[luigi1973]

if i have understood correctly, imap and map provide forward and inverse functionality on CPU. When does Vispy use them? Grazie

[djhoese]

You are correct, it is the CPU equivalent of the GLSL. It isn't used much in the python code from what I can tell, but I did use it a lot for debugging things in my own application. It is nice to have a way to easily debug transformation issues when you have a python equivalent. This is another reason I would like if we used PROJ.4 definitions for geographic transformations because the python code could just use pyproj.

[djhoese]

I'm removing this from the 0.7 milestone. While I still find this feature very important it will take a lot of work that we just don't have the time for to be included in 0.7. We've had renewed interest from @rougier and other original maintainers of VisPy which may influence how and where these definitions come from (glumpy? my SIFT project?) and who implements them. Our roadmap that I've put together also points out our need to get away from OpenGL 2.0 limitations and eventually OpenGL entirely. Not being limited to OpenGL 2.0 may mean it would be best to have these transforms not only occur in the vertex shader but also in a tessellation shader to get the "best" end result.

Additionally, there has been a lot of development in pyproj and the PROJ C library over the last couple years. PROJ is now very strict about CRS definitions and what is possible. We'll need to think about if we just want our transforms to be "good enough" or if we want to allow for something as complex as what PROJ has implemented (time-dependent CRS/datums, etc).

[kmuehlbauer]

@djhoese There is this interesting paper (https://doi.org/10.14311/gi.15.1.5) which compares GDAL warping with OpenCL based warping. It concludes that not in all cases OpenCL is the better approach. The paper is already 5 years old, but the bottom line should still hold. Not sure if reimplementation of geographical transformations in VisPy is the right choice for now. I'm all good with removing from the 0.7 milestone.

[djhoese]

@kmuehlbauer That is interesting. I didn't read beyond the abstract, but I think the major thing about my comment above is the coordinates that PROJ is able to produce and their accuracy. There has been some talk about OpenCL or CUDA versions of PROJ algorithms in their mailing list, but I'm not sure it really went anywhere.

[kmuehlbauer]

@djhoese Should we talk about this at the next dev meeting?

[djhoese]

Sure. If you can update the meeting notes and include that in the agenda that would be great (and please link to this issue).

[djhoese]

As mentioned in #820, we should try making an example that involves cartopy (if possible). Here is the glumpy example: https://github.com/glumpy/glumpy/blob/master/examples/cartopy.py

[djhoese]

In our developer meeting today we discussed this issue and decided that this will not be included as part of vispy. The pyproj (and underlying PROJ C library) are much too complicated to expect vispy to include an equivalent set of projections with equivalent precision. You run into issues whether you transform on the CPU or the GPU with only transforming some of the vertices instead of being able to interpolate along a line.

We recommend users transform coordinates with pyproj before giving them to a vispy visual or that you write custom transform objects with your own shaders for doing the transformation. If someone is really interested in attempting to port the C versions of these transformations to shader language (GLSL 1.20 compatible or not) and put them under the vispy "umbrella" then let me know and we can create a "vispy-geo" or similar project.

Now all of this said, I have an application where I created my own custom projection transformation shaders and I depend on the transformation being in the GPU (at least with how I'm doing it now). If I was to transition to not using these transformations in the GPU and instead did everything in the CPU, I'm not sure how it would effect my application. So my opinion may change about this topic in the future but for now we're saying it won't be supported. I will attempt to find time to make a simple example of using pyproj with some data and plot it on different projections, but I'm not sure when that will happen. I will leave this open until that happens though.

[djhoese]

I also wanted to add that @kmuehlbauer pointed out a paper about OpenCL (note: CL not GL) resampling of data versus gdalwarp CPU-based warping. These geographic transformations being discussed here are purely transforming the vertices and we depend on the GPU stretching the data between two vertices.

Another option in addition to what I said in my previous comment would be to have a pure OpenGL or CUDA implementation of these algorithms and share the results with VisPy (see #210 or #1800).

[jni]

@djhoese

I have an application where I created my own custom projection transformation shaders and I depend on the transformation being in the GPU

Mind adding a link to the code? 😂

We've had a few napari users request geographical transforms and we need to fiddle with our own transform machinery but having a few good examples on the VisPy side would be extremely useful. Currently we only have vismap as an example which only implements one very simple projection.

[djhoese]

This is very very ugly and was written when PROJ4 was more acceptable (it is generally frowned upon to depend on PROJ4 strings/dicts in favor of WKT):

https://github.com/ssec/sift/blob/master/uwsift/view/transform.py#L636