Ideas for additional checks for "scratches" or shapfile lines crossing the image

[djhoese]

Problem description

There are various cases where due to the nature of shapefiles in lon/lat space and images in various projections' space, that you get drawn lines from shapefiles that go straight across the image. This essentially comes from "point N" of a polygon mapping to the outside left side of the image and "point N +1" mapping to the outside right side of the image and pycoast drawing a line connecting these two points. This often happens because the points straddle the anti-meridian of the target projection or the anti-meridian of the lon/lat space (-180/180) depending on the shapefiles used.

@lobsiger had done some work in the past to minimize the number of these lines in #59 and I think a couple other PRs. I'm starting this issue to talk about other options that weren't implemented yet and I think could easily be added, but may affect performance.

The Code

The main part of the code that would be modified is this method here:

pycoast/pycoast/cw_base.py

Line 1088 in 4a34299

if (x >= 1e30).any() or (y >= 1e30).any():

This line is essentially a boolean mask used to say "is this point valid in this projection" and is used in a later line to determine where to split a polygon into smaller segments to avoid drawing a line from point N to point N+1 where it would draw a large non-real line across the image.

Adding any additional conditions to this should be easy as we can just do a boolean AND or OR with a numpy array.

Possible solutions

My first thought was to take the difference between each coordinate point in the X dimension (the only dimension that has wrapping/non-contiguous coordinates) using np.diff and look for anything larger than the X-extents of the image. I realized later that this could be a problem for very small geographic regions where the points for larger features like country borders and coastlines might be much further apart than looking at a single province or city, even though you'd still want the line to show up.

My next thought was still to do the np.diff but to maybe look for any difference larger than X degrees/meters (depending on the output projection). A good value for this could be determined offline (not at runtime, just store it in a constant) by taking all the GSHHG shapefiles/polygons and determine the largest difference between their points in geocentric meters. This should give a clear indication of "point N is at location A, but point N+1 is at location A+10000000 meters" and work for global images and local (city-level) images. @lobsiger's existing checks should limit polygons that are completely outside the image already, but this will handle the anti-meridian case...I hope.

Thoughts @mraspaud @lobsiger ?

[mraspaud]

Could we split the polygons somehow with the antimeridian line?

[djhoese]

Well the shapefiles we use currently are split at the -180/180 longitude anti-meridian. To split them at the output image's projection's anti-meridian (if one exists) I think we'd do a similar computation to what I'm talking about, but we'd have to come up with a consistent way of determining the circumference of the projection. In the past I've cheated and done 2x the X distance from (0E, 0N) and (180E, 0N). I'm not sure that is reliable for all projections.

So I think my original hope for a solution still applies, but I think I can get rid of the geocentric meters idea. We find the largest difference between two points in the shapefiles in degrees (offline, one time, not at run time). We then, at runtime, transform that number of degrees at our image's origin and take the max of the x/y distance. If a polygon we're about to draw has two vertices that are further than that then we can assume the polygon is connecting across a "jump" in the projection (anti-meridian or other non-contiguous portion).

[lobsiger]

A rubber 3D sphere must be cut somehow to be distorted/flattened/projected to 2D. As I have mentioned before depending on the type of projection at hand there are cuts along meridians (direct), along the equator (transverse) or along other special great circles (oblique). Exotic double periodic projections usually use cross like cuts. The conformal and widely used 'stere' azimuthal projection just uses an opposite anti pinhole enlarged to infinity. Other azimuthal projections enlarge the anti pinhole to a limiting circle. I doubt that there is a solution that handles all these rather different cases in an elegant and fast way.

Proj has even much more exotic projections but looking at the satpy distributed 'areas.yaml' there are in fact surprisingly few of them in use. Some potentially useful pseudo conic Proj projections cannot be handled with the current resamplers we have.

PyCoast has a fast way of excluding shapes (of all existing GSHHG resolutions) that might be harmful. The idea is to simply avoid crossing of the above mentioned cuts which also leaves the possibility to fill closed shapes intact. Most of the remaining problems fall in one of these two categories:

A)
Limited areas of direct aspect cylindrical or conical projections where the huge Eurasia polygon can be a problem. There is usually a simple workaround by shifting lon_0 and using asymmetric area_extents. There is no corresponding problem with transverse aspects (e.g. for all sorts of UTM) as the Eurasia polygon does not cross the equator. There might be very rare cases (I have seen that at least twice) where images from LEOs using 'omerc_bb' will have a problem with the Eurasia shape.

B)
Full area of the globe where users insist to use cylindrical, conic or pseudo conic projections with anti meridians different from +/-180° or even ask for very special projection aspects by adding an 'ob_tran'.

https://proj.org/en/9.2/operations/projections/ob_tran.html

It must be noted that plotting the grid will also fail in PyCoast. Furthermore a clean solution would not just lift the pen when crossing one of these cuts but introduce two points as close as possible on both sides of the cut to draw the coastlines (and also the lat/lon grid lines) to the 'circumference' of the projection.

While it is always possible to choose deliberately bad parameters to demonstrate the "problem" (even for azimuthal projections) I see NO common use case that cannot easily be fixed for A) and just recommend to let lon_0 = 0° for B) as most users will do anyway. Bottom line: "If it ain't broken don't fix it!"

Ernst

[djhoese]

and just recommend to let lon_0 = 0°

Asking users to change the projection of their data/image is not a solution in my opinion. And while we may not find a solution that works for every projection I think we can limit the number times we see issues.

I think the solution I'm hoping to implement regarding this issue is one for non-global images where lines are going across the image and would otherwise not be in the image. I think this in addition to the existing filtering and segmenting behavior should stop huge lines across the image, but could result in global images having polygons on the edge of the image/projection that just stop short of the edge (as we break the polygon at that segment).

[lobsiger]

I think the solution I'm hoping to implement regarding this issue is one for non-global images where lines are going across the image and would otherwise not be in the image. @djhoese can you show us some common satellite image use cases that still have that problem (and cannot easily be fixed with adjusting lon_0)? Have satpy users still reported such scratch problems after my #59 fix?

[djhoese]

Adjusting lon_0 is not a solution. It requires me to change the projection that I have and the way my image looks. The projection being used or where data is located may not be something a user can change.

We already talked about cases where this is an issue in #98. In that case it was a prime meridian shift on the lon/lat projection.

[lobsiger]

@djhoese my Honolulu picture in #98 shows the Eurasia problem that has the simple workaround. As explained you change the prime meridian lon_0=170° and shift the area extent by 10°. Then using the same projection your picture looks exactly the same with the same data used but without scratches. You already said you do not like any workaround but want a general solution. My point is that users editing their own areas.yaml should be aware of the Eurasia problem and generally know what they (can) do.

[djhoese]

Some users can't change their projection. It is out of their control. Changing the projection is not a solution. I suppose these are always going to generate PNGs so there is no CRS/WKT information so maybe it doesn't matter. Bottom line is users shouldn't have to modify their projection/area.

[lobsiger]

If pycoast users can't do anything they will have some admin that I hope knows what he does. This guy will not propose/prepare projections that produce scratches all over the place. I agree that we face an interesting problem to solve but can't see that it is an important problem pycoast has. It's rather a feature than a bug. In any case to me it is of far less importance than the resampling problems/bugs also mentioned in #98. And yes, my images are just *.png, *.jpg and *.webp that contain no geo reference info.

[djhoese]

If pycoast users can't do anything they will have some admin that I hope knows what he does. This guy will not propose/prepare projections that produce scratches all over the place.

Everyone has different requirements. Your use case is only your use case. Some people are using area definitions that are configured by someone else, or used by some other tool, or are converted from some other geographic definition. Some users have to give their images/geotiffs to another third-party tool that will assume certain things. Some people don't have the control or flexibility with changing things that you have. Please accept that you are not the only user and your use case is not the only use case. If someone was using pycoast and saw that the image for data case 1 produces perfectly fine coastlines and for data case 2 it produces "scratches" over the image then it should be clear to them that this is a bug in the tool they are using (pycoast).

As far as resampling problems, I mentioned it in the mailing list I think, but we need to put together some github issues on the pyresample repository if they don't already exist. Some of the tests were using bilinear and gradient_search which I'm not sure are very robust. Gradient has had some updates in the past 6 months or so, so it may behave a little better, but I also wouldn't be surprised if it didn't. Some of the issues were specific to the thin version of the MODIS data if I recall correctly, but yeah definitely something we need to track down.

[lobsiger]

@djhoese I do not say that my use case is something special and of course I cannot know all the PyCoast use cases. I just say that as a SatPy user I'm quite happy with PyCoast as is and have not seen other users complaining about the scratch problem lately.

Two things we probably agree is that:

1) A better handling of coastline shape files should not slow down PyCoast considerably.
2) The possibility to fill islands and continents should be preserved. This comes in handy
e.g. in Ocean Science and I have used it for ocean color or temperature data myself.

If the prime meridian is chosen in the middle of e.g. Mercator images the huge Eurasia polygon can not only scratch from left to right but also impair filling. My easy workaround fixes that as well (example attached). If a solution for the PM in the center only lifts the pen, then the filling problem might remain. If @mraspaud this solution actually cuts the polygon in two or more polylines then filling is not possible anymore unless the fragments are again closed -- eventually inserting even multiple points -- along the projection "cut" (the anti meridian in below Mercator case used in NOAA Ocean Prediction Center MSLP charts) on each sides.