{"action":"created","issue":{"url":"https://api.github.com/repos/spacetelescope/spherical_geometry/issues/125","repository_url":"https://api.github.com/repos/spacetelescope/spherical_geometry","labels_url":"https://api.github.com/repos/spacetelescope/spherical_geometry/issues/125/labels{/name}","comments_url":"https://api.github.com/repos/spacetelescope/spherical_geometry/issues/125/comments","events_url":"https://api.github.com/repos/spacetelescope/spherical_geometry/issues/125/events","html_url":"https://github.com/spacetelescope/spherical_geometry/issues/125","id":274708827,"number":125,"title":"Bugs in area and intersection","user":{"login":"mcdigman","id":1415033,"avatar_url":"https://avatars0.githubusercontent.com/u/1415033?v=4","gravatar_id":"","url":"https://api.github.com/users/mcdigman","html_url":"https://github.com/mcdigman","followers_url":"https://api.github.com/users/mcdigman/followers","following_url":"https://api.github.com/users/mcdigman/following{/other_user}","gists_url":"https://api.github.com/users/mcdigman/gists{/gist_id}","starred_url":"https://api.github.com/users/mcdigman/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/mcdigman/subscriptions","organizations_url":"https://api.github.com/users/mcdigman/orgs","repos_url":"https://api.github.com/users/mcdigman/repos","events_url":"https://api.github.com/users/mcdigman/events{/privacy}","received_events_url":"https://api.github.com/users/mcdigman/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"assignees":[],"milestone":null,"comments":1,"created_at":"2017-11-17T00:11:05Z","updated_at":"2018-02-07T21:06:05Z","closed_at":null,"author_association":"NONE","body":"The following code demonstrates two bugs which may or may not be related. The first is that area() for a valid polygon gives a negative answer, which is not possible. The second bug is that the intersection of the two geometries in the test case gives the complement of the intersection instead. Both bugs exists in v1.0.12, and also in a clone from the github made this afternoon. In this test case, whether the bugs are present or not appears to depend unpredictably on minor changes in the geometries which geometries it is present for seems to depend on the version. I am using Python 2.7.14, Matplotlib 2.1.0, numpy 1.13.3, scipy 1.0.0, and astropy 2.0.2 on OSX 10.11.6.\r\n\r\n\r\n`\r\nimport numpy as np\r\nfrom spherical_geometry.polygon import SphericalPolygon\r\nimport spherical_geometry.vector as sgv\r\n\r\nn_fill = 10\r\n\r\ntheta1_high_fill = np.full(n_fill,5.)\r\ntheta1_low_fill =np.full(n_fill, -65.)\r\nphi1_high_fill = np.linspace(50.-360.,50.-1.,n_fill)\r\nphi1_low_fill = phi1_high_fill[::-1]\r\ntheta1s = np.hstack([theta1_high_fill,theta1_low_fill,theta1_high_fill[0]])\r\nphi1s = np.hstack([phi1_high_fill,phi1_low_fill,phi1_high_fill[0]])\r\n\r\nphi_in1 = 0.\r\ntheta_in1 = 0.\r\n\r\nbounding_xyz1 = np.array(sgv.radec_to_vector(phi1s,theta1s,degrees=True)).T\r\ninside_xyz1 = np.array(sgv.radec_to_vector(phi_in1,theta_in1,degrees=True))\r\nsp_poly1 = SphericalPolygon(bounding_xyz1,inside=inside_xyz1)\r\n\r\ntheta2_high_fill = np.full(n_fill,-30.)\r\ntheta2_low_fill =np.full(n_fill, -85.)\r\nphi2_high_fill = np.linspace(65-360.,55.-10.,n_fill)\r\nphi2_low_fill = np.linspace(5-360.,5.-10.,n_fill)[::-1]\r\ntheta2s = np.hstack([theta2_high_fill,theta2_low_fill,theta2_high_fill[0]])\r\nphi2s = np.hstack([phi2_high_fill,phi2_low_fill,phi2_high_fill[0]])\r\n\r\nphi_in2 = 0.\r\ntheta_in2 = -70.\r\n\r\nbounding_xyz2 = np.array(sgv.radec_to_vector(phi2s,theta2s,degrees=True)).T\r\ninside_xyz2 = np.array(sgv.radec_to_vector(phi_in2,theta_in2,degrees=True))\r\nsp_poly2 = SphericalPolygon(bounding_xyz2,inside=inside_xyz2)\r\n\r\nint_poly = sp_poly2.intersection(sp_poly1)\r\n\r\ntheta_out = 30.*np.pi/180.\r\nphi_out = 0.\r\noutside_xyz = np.array(sgv.radec_to_vector(phi_out,theta_out,degrees=False))\r\n\r\ntheta_in3 = -40.*np.pi/180.\r\nphi_in3 = 0.\r\ninside_xyz3 = np.array(sgv.radec_to_vector(phi_in3,theta_in3,degrees=False))\r\n\r\nprint \"area int, area 1,area 2: \",int_poly.area(),sp_poly1.area(),sp_poly2.area()\r\n\r\nprint \"Should be False: \",int_poly.contains_point(outside_xyz)\r\n\r\n\r\nfrom mpl_toolkits.basemap import Basemap\r\nimport matplotlib.pyplot as plt\r\nm = Basemap(projection='moll',lon_0=0)\r\n\r\nfig = plt.figure(figsize=(10,5))\r\nax = fig.add_subplot(111)\r\n\r\nsp_poly1.draw(m,color='red')\r\nsp_poly2.draw(m,color='blue')\r\nint_poly.draw(m,color='green')\r\nplt.show()`\r\n\r\nThe following code demonstrates that the intersection is in fact the complement:\r\n`print \"Should be False,False,False: \",int_poly.contains_point(outside_xyz),sp_poly1.contains_point(outside_xyz),sp_poly2.contains_point(outside_xyz)\r\nprint \"Should be False,True ,False: \",int_poly.contains_point(inside_xyz1),sp_poly1.contains_point(inside_xyz1),sp_poly2.contains_point(inside_xyz1)\r\nprint \"Should be False,False,True : \",int_poly.contains_point(inside_xyz2),sp_poly1.contains_point(inside_xyz2),sp_poly2.contains_point(inside_xyz2)\r\nprint \"Should be True ,True ,True : \",int_poly.contains_point(inside_xyz3),sp_poly1.contains_point(inside_xyz3),sp_poly2.contains_point(inside_xyz3)\r\n`\r\n\r\nThe attached picture shows the 3 geometries in question, and that the 'inside' point for the intersection is in fact not inside the intersection region.\r\n![sph_geo_bug1](https://user-images.githubusercontent.com/1415033/32922609-031bd402-cb01-11e7-984a-8a7ea58f8a0b.png)\r\n"},"comment":{"url":"https://api.github.com/repos/spacetelescope/spherical_geometry/issues/comments/363910489","html_url":"https://github.com/spacetelescope/spherical_geometry/issues/125#issuecomment-363910489","issue_url":"https://api.github.com/repos/spacetelescope/spherical_geometry/issues/125","id":363910489,"user":{"login":"mcdigman","id":1415033,"avatar_url":"https://avatars0.githubusercontent.com/u/1415033?v=4","gravatar_id":"","url":"https://api.github.com/users/mcdigman","html_url":"https://github.com/mcdigman","followers_url":"https://api.github.com/users/mcdigman/followers","following_url":"https://api.github.com/users/mcdigman/following{/other_user}","gists_url":"https://api.github.com/users/mcdigman/gists{/gist_id}","starred_url":"https://api.github.com/users/mcdigman/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/mcdigman/subscriptions","organizations_url":"https://api.github.com/users/mcdigman/orgs","repos_url":"https://api.github.com/users/mcdigman/repos","events_url":"https://api.github.com/users/mcdigman/events{/privacy}","received_events_url":"https://api.github.com/users/mcdigman/received_events","type":"User","site_admin":false},"created_at":"2018-02-07T21:06:05Z","updated_at":"2018-02-07T21:06:05Z","author_association":"NONE","body":"I pulled the most recent version of the package this afternoon. I see that my original test case has been fixed. However if I change nothing about the geometry except moving all the theta_in points to the northern hemisphere, the code still misbehaves:\r\n```\r\nimport numpy as np\r\nfrom spherical_geometry.polygon import SphericalPolygon\r\nimport spherical_geometry.vector as sgv\r\n\r\nn_fill = 10\r\n\r\ntheta1_high_fill = np.full(n_fill,5.)\r\ntheta1_low_fill =np.full(n_fill, -65.)\r\nphi1_high_fill = np.linspace(50.-360.,50.-1.,n_fill)\r\nphi1_low_fill = phi1_high_fill[::-1]\r\ntheta1s = np.hstack([theta1_high_fill,theta1_low_fill,theta1_high_fill[0]])\r\nphi1s = np.hstack([phi1_high_fill,phi1_low_fill,phi1_high_fill[0]])\r\n\r\nphi_in1 = 0.\r\ntheta_in1 = 70.\r\n\r\nbounding_xyz1 = np.array(sgv.radec_to_vector(phi1s,theta1s,degrees=True)).T\r\ninside_xyz1 = np.array(sgv.radec_to_vector(phi_in1,theta_in1,degrees=True))\r\nsp_poly1 = SphericalPolygon(bounding_xyz1,inside=inside_xyz1)\r\n\r\ntheta2_high_fill = np.full(n_fill,-30.)\r\ntheta2_low_fill =np.full(n_fill, -85.)\r\nphi2_high_fill = np.linspace(65-360.,55.-10.,n_fill)\r\nphi2_low_fill = np.linspace(5-360.,5.-10.,n_fill)[::-1]\r\ntheta2s = np.hstack([theta2_high_fill,theta2_low_fill,theta2_high_fill[0]])\r\nphi2s = np.hstack([phi2_high_fill,phi2_low_fill,phi2_high_fill[0]])\r\n\r\nphi_in2 = 0.\r\ntheta_in2 = 70.\r\n\r\nbounding_xyz2 = np.array(sgv.radec_to_vector(phi2s,theta2s,degrees=True)).T\r\ninside_xyz2 = np.array(sgv.radec_to_vector(phi_in2,theta_in2,degrees=True))\r\nsp_poly2 = SphericalPolygon(bounding_xyz2,inside=inside_xyz2)\r\n\r\nint_poly = sp_poly2.intersection(sp_poly1)\r\n\r\ntheta_out = -50.*np.pi/180.\r\nphi_out = 0.\r\noutside_xyz = np.array(sgv.radec_to_vector(phi_out,theta_out,degrees=False))\r\n\r\ntheta_in3 = 70.*np.pi/180.\r\nphi_in3 = 0.\r\ninside_xyz3 = np.array(sgv.radec_to_vector(phi_in3,theta_in3,degrees=False))\r\n\r\nprint \"area int, area 1,area 2: \",int_poly.area(),sp_poly1.area(),sp_poly2.area()\r\nprint \"Should be True , True :\",int_poly.area()<=sp_poly1.area(),int_poly.area<=sp_poly2.area()\r\n\r\nprint \"Should be False,False,False: \",int_poly.contains_point(outside_xyz),sp_poly1.contains_point(outside_xyz),sp_poly2.contains_point(outside_xyz)\r\nprint \"Should be True ,True ,True : \",int_poly.contains_point(inside_xyz1),sp_poly1.contains_point(inside_xyz1),sp_poly2.contains_point(inside_xyz1)\r\nprint \"Should be True ,True ,True : \",int_poly.contains_point(inside_xyz2),sp_poly1.contains_point(inside_xyz2),sp_poly2.contains_point(inside_xyz2)\r\nprint \"Should be True ,True ,True : \",int_poly.contains_point(inside_xyz3),sp_poly1.contains_point(inside_xyz3),sp_poly2.contains_point(inside_xyz3)\r\n\r\ninside_res = list(int_poly.inside)\r\nfor itr in xrange(0,len(inside_res)):\r\n    print \"Should be True ,True ,True : \",int_poly.contains_point(inside_res[itr]),sp_poly1.contains_point(inside_res[itr]),sp_poly2.contains_point(inside_res[itr])\r\n\r\n\r\nfrom mpl_toolkits.basemap import Basemap\r\nimport matplotlib.pyplot as plt\r\nm = Basemap(projection='moll',lon_0=0)\r\n\r\nsp_poly1.draw(m,color='red')\r\nsp_poly2.draw(m,color='blue')\r\nint_poly.draw(m,color='green')\r\nplt.show()\r\n```\r\nGives:\r\n```\r\narea int, area 1,area 2:  6.7868081724 6.27810880223 9.69856292105\r\nShould be True , True : False False\r\nShould be False,False,False:  True False False\r\nShould be True ,True ,True :  False True True\r\nShould be True ,True ,True :  False True True\r\nShould be True ,True ,True :  False True True\r\nShould be True ,True ,True :  True False False\r\n```\r\nAnd produces the attached plot. The intersection routine is choosing an inside point which is not inside either of the original geometries, and as far as I can tell, int_poly is the complement of the correct intersection. \r\n\r\n![sph_geo_bug2](https://user-images.githubusercontent.com/1415033/35940735-208c8e7e-0c1e-11e8-9a16-ba10c2d4f7d9.png)\r\n"}}
{"actor":{"display_login":"mcdigman"}}