OSError: broken data stream when reading image file

the remedy is to manually change the handler to be able to read truncated files; inside the file pydicom/pixel_data_handlers/pillow_handler.py explicitly set ImageFile to be able to read truncated images, i.e. change from:

try:
    import PIL
    from PIL import Image, features

to

try:
    import PIL
    from PIL import Image, features, ImageFile
    ImageFile.LOAD_TRUNCATED_IMAGES = True

then you’ll get to read the pixel array no matter what, but you’re still fucked if your actual image is corrupt like this:

CuPy has a builtin experimental profiler (cupyx.time) that can accurately assess the GPU runtime. But for simplicity the benchmark here only measures the elapsed time, and represented by the mean value based on multiple iterations. The overhead of the first invocation is excluded but the transfer to the host is not.

While this is all impressive, it is worth noting that NumPy inversions can be made faster (a factor of few) by skipping identity matrix recreation (directly call np.linalg.solve) or by directly calling the barebone lapack routine.

All calculations were performed on a Google Colab instance with:

CPU: Intel(R) Xeon(R) @ 2.00GHz
GPU_0: Tesla T4

numpy == 1.19.5
cupy == 9.1.0
cuda == 11.0.221

with x and y being the coordinates of n vertices. If you have the coordinates as one-dimensional arrays x and y, using the numpy.roll() function this can be implemented as a beautiful one-liner:

np.abs(np.dot(x, np.roll(y, 1)) - np.dot(y, np.roll(x, 1))) / 2.

The plain map depicting just land and water is obtained through the ‘naturalearth_lowres’ dataset (don’t think I saw one for coastline). It returns boundaries for each country, but it can be suppressed by assigning the same color to both color and edgecolor for each polygon. Since it would be easiest to get the final fraction by counting pixel intensities in a 2/3D image array, we would need to produce the image arrays based on an equal-area projection (like Mollweide).

While creating different projections is easy in geopandas, I couldn’t find a reliable way of shifting the center of a projection only except in Mercator. And so my work-around was to shift the projection first in Mercator for a given angle, and then re-project in Mollweide. Here I note that when assembling a dataframe from scratch its coordinate reference system (crs) has to be explicitly assigned, especially if you later wish to re-project in a different system (Mercator=’EPSG:4326’, Mollweide=’ESRI:54009’).

Also I couldn’t find a way to directly export plotted geopandas objects into a 2/3D array. My workaround was to literally save the plot without any padding and re-read the image as numpy array. We lose the native coordinates, but it doesn’t matter for our case since we will be just counting pixels by their intensities.

Overall, the script creates two Mollweide projections where one of them is shifted by 180 degrees and co-axially flipped. Each projection is ultimately represented as a grayscale 2D array with land and water having uniform intensities of 0.5 and 0 respectively. When these are added together, the combined image will have pixel values of 1 wherever antipodal points both occur on land. Running the script gives:

probability of antipode being on land: 0.150

and according to wiki it is indeed 15%; in short, it is the consequence of most of the land mass being clumped in one hemisphere.

MSE = 0.00010