warnings of invalid value encountered in true_divide and invalid value encountered in double_scalars in

[DikraK]

Hi all,

My objective is to found if a satellite scene (sentinel-1) intersects with a region defined by a polygon in a shapefile.

Here is my script:

#!/usr/bin/env python

import numpy as np
from satpy.utils import debug_on
from pyresample.spherical import SphPolygon
import shapefile
from zipfile import ZipFile
import xml.etree.ElementTree as ET

def get_sphpolygon(shp_file):
    """Generate  pyresample SphPolygon from a shapefile"""
    sf           = shapefile.Reader(shp_file)
    feature  = sf.shape(0).__geo_interface__
    vertices = feature["coordinates"][0]

    return SphPolygon(np.deg2rad(vertices))

def get_swath_sphpolygon(s1_file):
    """return S1 swath polygon"""

    try:
        with ZipFile(s1_file, "r") as zipobj:
            for name in zipobj.namelist():
                if name.endswith("preview/map-overlay.kml"):
                    root            = ET.fromstring(zipobj.read(name))
                    coordinates     = root.findall(".//coordinates")
                    if len(coordinates) > 0:
                        latlon_quad = coordinates[0].text
                        points_lonlat = [
                            (float(latlon.split(",")[0]), float(latlon.split(",")[1]))
                            for latlon in latlon_quad.split()
                        ]

                        scan_polygon = SphPolygon(np.deg2rad(points_lonlat))

                        return scan_polygon
    except BadZipFile:
        logger.info("Found bad zip file {}, skip processing".format(s1_file))
        pass

    return None


def is_intersected(s1_file, polygon_region):
    """check if the swath is intersected with the target region"""
    scan_polygon = get_swath_sphpolygon(s1_file)

    return (
        scan_polygon
        and scan_polygon._is_inside(polygon_region)
        or scan_polygon.intersection(polygon_region)
    )


## Apply the above functions to found if there are intersections

# 1. transform the shapefile polygon  to SphPolygon
region_polygons = get_sphpolygon(shp_file_name)  # shp_file_name : is the name of the shapefile

# 2. check if the sentinel-1 scenes intersect with the polygon
swath = 0 
for f in files: # is a list of sentinel 1 filename
       try:
            with ZipFile(f):
                if is_intersected(f, region_polygons):
                    swath += 1         
        except BadZipFile as e:
            print(e)
            pass
    else:
        pass

My script run just fine but I always have a warning at the beginning of the run:

.conda/envs/myenv/lib/python3.8/site-packages/pyresample/spherical.py:119: RuntimeWarning: invalid value encountered in true_divide
  self.cart /= np.sqrt(np.einsum('...i, ...i', self.cart, self.cart))
.conda/envs/myenv/lib/python3.8/site-packages/pyresample/spherical.py:171: RuntimeWarning: invalid value encountered in double_scalars
  return (val + mod) % (2 * mod) - mod

It seems to be caused by the lines:

 return (
        scan_polygon
        and scan_polygon._is_inside(polygon_region)
        or scan_polygon.intersection(polygon_region)
    )

Any hints on why this warning keeps showing up ?

Environment Info:

• OS: ubuntu 18
• Satpy Version: 0.29.0
• PyResample Version: 1.21.0
• Python Version : 3.8.6
• Others: working in a conda environment

[djhoese]

Sorry for the late response on this, I thought I had answered when you first made it. These types of warings are typical when NaNs are involved. For normal operations where you're dealing with an entire swath this isn't too big of a deal, but if they are showing up in your polygon or a bounding box definition then that probably means things aren't being described properly. Would it be possible for you to print out points_lonlat and paste the output here?

Edit: Also I should point out that you import satpy, but it doesn't look like you ever use it. I guess this is more of a pyresample/numpy question than satpy. No big deal though.

[DikraK]

That's ok! thanks for your response.

if I print out the points_lonlat for a given Sentinel-1 image, I have:

points_lonlat : [(-65.055031, 58.893234), (-57.768929, 59.644321), (-56.693855, 55.808479), (-63.270271, 55.089134)]

And here is the polygon_region printed out:

[[-113.21558887   48.70120642]
 [-115.66462668   52.3475516 ]
 [-118.01570297   53.51220513]
 [-117.66317189   53.85272687]
 [-116.78904762   54.72342275]
 [-113.33540294   55.39778732]
 [-106.41200236   55.3383318 ]
 [-104.29057378   55.69794751]
 [-102.16907227   55.91592497]
 [ -99.42814784   56.12728191]
 [ -97.16086887   56.576744  ]
 [ -94.95566517   57.02997474]
 [ -93.51207202   57.58351759]
 [ -92.23912736   57.28988879]
 [ -91.81481247   57.02690295]
 [ -92.45277624   56.64276048]
 [ -93.1080232    56.30259784]
 [ -94.37275656   55.83148566]
 [ -95.37462313   55.16534632]
 [ -96.43904132   54.45761652]
 [ -95.28331232   53.1090265 ]
 [ -92.60817156   52.43657916]
 [ -90.17816932   50.85042486]
 [ -89.17799851   48.43713542]
 [ -92.23912736   46.92733286]
 [ -94.20916077   46.80299942]
 [ -95.08809875   45.20315943]
 [ -97.78552912   45.20315943]
 [ -99.05847379   47.19573433]
 [ -99.66463791   47.03072411]
 [-102.81669138   48.47733494]
 [-113.21558887   48.70120642]]

It is really the intersection and _is_inside that raise the warning.

[djhoese]

Ok so you're data is fine. 😕

It has been a little while since I've looked at the SphPolygon code, but I think some of the calculations end when values become invalid so maybe this is expected 🤷‍♂️ It is possible to tell numpy to not raise these and catch any regular python RuntimeWarnings (which I think is what these are) by doing:

import warnings
import numpy as np

with np.errstate(invalid="ignore"), warnings.catch_warnings():
    # your code that raises the warnings

It isn't great, but it may be necessary in this case. @mraspaud any thoughts on this?