Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Location does not return a country code when truncated? #12

Closed
Zamirk opened this issue Jan 13, 2023 · 6 comments
Closed

Location does not return a country code when truncated? #12

Zamirk opened this issue Jan 13, 2023 · 6 comments

Comments

@Zamirk
Copy link

Zamirk commented Jan 13, 2023

The following Latitude and Longitude (45.8, 16) is within the borders of Croatia. (HR)

For some reason, it does not return a country code. I think it occurs for all locations with (16) Longitude.

Interestingly, adding 0.00001 to the 16 value fixes the issue.

All 3 datasets return the same results

  • boundaries180x90.ser
  • boundaries60x30.ser
  • boundaries360x180.ser

Scala Sample. I used the method call

boundaries.getIds(16, 45.8)

for {
         fileStream <- F.delay(new FileInputStream(config.serFile))
         boundaries <- F.delay(CountryBoundaries.load(fileStream))
         result     <- F.delay(boundaries.getIds(longitude, latitude).asScala.toList)
         _          <- F.delay(fileStream.close())
} yield result
@westnordost
Copy link
Owner

westnordost commented Jan 13, 2023

Which version of the lib are you using? Is this also reproducible for latitude? What about longitude 17 (etc)?

So this might have something to do with that at the edges of the raster cells, a position would not be considered to be inside any of the two neighboring raster cells. However, lon 16 is not even at a raster edge for 2 of the 3 data sets (360 / 60 = 6, so the edges are at lon 0, 6, 12, 18,...)

@westnordost
Copy link
Owner

Can reproduce. This is how the cell in which the point in polygon check is done looks like:

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [15.74585,45.0638], [15.78842,45.11519], [15.76371,45.16508], [15.83512,45.22459], [15.98412,45.23088], [16.0,45.2153109], [16.0,45.0], [15.7741511,45.0], [15.74585,45.0638]
          ]
        ]
      },
      "properties": {
        "country": "BA"
      }
    },
    {
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [15.0,45.0], [15.0,45.4953588], [15.02385,45.48533], [15.05187,45.49079], [15.16862,45.42309], [15.27758,45.46678], [15.33051,45.45258], [15.38188,45.48752], [15.30249,45.53224], [15.29837,45.5841], [15.27747,45.60504], [15.31027,45.6303], [15.34695,45.63382], [15.34214,45.64702], [15.38952,45.63682], [15.4057,45.64727], [15.34919,45.71623], [15.30872,45.69014], [15.25423,45.72275], [15.40836,45.79491], [15.47531,45.79802], [15.47325,45.8253], [15.52234,45.82195], [15.57952,45.84953], [15.64185,45.82915], [15.66662,45.84085], [15.70411,45.8465], [15.68232,45.86819], [15.68383,45.88867], [15.67967,45.90455], [15.70636,45.92116], [15.7032759,46.0], [16.0,46.0], [16.0,45.2153109], [15.98412,45.23088], [15.83512,45.22459], [15.76371,45.16508], [15.78842,45.11519], [15.74585,45.0638], [15.7741511,45.0], [15.0,45.0]
          ]
        ]
      },
      "properties": {
        "country": "HR"
      }
    },
    {
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [15.0,45.4953588], [15.0,46.0], [15.7032759,46.0], [15.70636,45.92116], [15.67967,45.90455], [15.68383,45.88867], [15.68232,45.86819], [15.70411,45.8465], [15.66662,45.84085], [15.64185,45.82915], [15.57952,45.84953], [15.52234,45.82195], [15.47325,45.8253], [15.47531,45.79802], [15.40836,45.79491], [15.25423,45.72275], [15.30872,45.69014], [15.34919,45.71623], [15.4057,45.64727], [15.38952,45.63682], [15.34214,45.64702], [15.34695,45.63382], [15.31027,45.6303], [15.27747,45.60504], [15.29837,45.5841], [15.30249,45.53224], [15.38188,45.48752], [15.33051,45.45258], [15.27758,45.46678], [15.16862,45.42309], [15.05187,45.49079], [15.02385,45.48533], [15.0,45.4953588]
          ]
        ]
      },
      "properties": {
        "country": "SI"
      }
    },
    {
      "type": "Feature",
      "geometry": {
        "type": "Point", 
        "coordinates": [16.0, 45.8]
      },
      "properties": {
        
      }
    }
    ]
}
Loading

It is directly on the right edge of the Croatia polygon in this cell.

Now, the winding number point in polygon algorithm used in this library considers only the lower and left edges of a polygon as inside of the polygon.

So, the actual bug is that this cell gets selected in the first place. The next cell to the right should be selected instead.

The algorithm used to get from a longitude to a raster cell X is, which seems correct:

Math.floor((180 + normalize(longitude, -180, 360)) / 360.0 * rasterWidth)

With rasterWidth = 360 and longitude = 16 this is

   Math.floor((180 + normalize(16, -180, 360)) / 360.0 * 360)
=> Math.floor(196 / 360.0 * 360)
=> Math.floor(196)
=> 196

However, due to floating point imprecision, the calculation inside the Math.floor ends up being 195.99999999999997. I believe this is the bug. Will investigate how to circumvent this.

@Zamirk
Copy link
Author

Zamirk commented Jan 14, 2023

Thank you for the prompt response and detailed explanation.

@westnordost
Copy link
Owner

So, I avoided the floating point imprecision and the resulting error by changing

 Math.floor((180 + normalize(longitude, -180, 360)) / 360.0 * rasterWidth) 

to

 Math.floor(rasterWidth * (180 + normalize(longitude, -180, 360)) / 360.0) 

(first multiply, then divide)


There are actually tests that are supposed to test the behavior exactly on those edges, but the tests used a rasterWidth of 2, so the issue was not triggered (because 0.5 * 2 is 1 even with floating point math).

@westnordost
Copy link
Owner

Will release an update now

@westnordost
Copy link
Owner

Done. Will take some time until it is available on mvn

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants