1

I am an absolute beginner with GIS data. I am attempting to write a tool to determine if a lat/long point resides within a polygon within a shapefile provided by SABS (https://nces.ed.gov/programs/sabs/). Basically, I am trying to write something that will tell me the primary school for a given address when internet is not available (no Google Maps, GreatSchools, etc.).

Anyway, I am using GeoTools with modified PointInPolygon code from https://github.com/ianturton/geotools-cookbook/blob/master/modules/filters/src/main/java/org/ianturton/cookbook/filters/PointInPolygon.java to search for lat/long points in polygons.

With that, I can see that there are lat/long points within the GeoTools suggested datafile (http://www.naturalearthdata.com/downloads/50m-cultural-vectors/). However, the SABS data does not report any points in polygons.

When I look at the polygon data for the two files, the points for the suggested data file are in the lat/long ranges I expect (+/-90, +/-180). The point data in the SABS file has values many magnitudes larger.

My guess is that the SABS data is in a projection that doesn't play well with querying if a point is inside the polygon. I don't know how to transform the data so that I can perform lat/long queries.

Can someone explain what I need to do to be able to perform a lat/long search with the SABS data?

SABS prj: PROJCS["USA_Contiguous_Albers_Equal_Area_Conic_USGS_version",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Albers"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-96.0],PARAMETER["Standard_Parallel_1",29.5],PARAMETER["Standard_Parallel_2",45.5],PARAMETER["Latitude_Of_Origin",23.0],UNIT["Meter",1.0]]

Natural Earth prj:

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

Code added to PointInPolygon to test data files:

for (double lat = -90.; lat < 90.; lat += 0.01) {
        for (double lon = -180; lon < 180; lon += 0.01) {
            final Point p = fac.createPoint(new Coordinate(lat, lon));
            final boolean flag = tester.isInShape(p);
            if (flag) {
                System.out.println(p + " is inside");
            }
        }
    }
  • you need to reproject one layer to match the other - see docs.geotools.org/stable/userguide/tutorial/geometry/… – Ian Turton Jan 24 '17 at 12:22
  • Thanks, @iant. I ran the CRSLab on the page to which you linked and was informed that the SABS data has no valid geometry, though I don't really know what to do with that information. – Todd Jan 24 '17 at 14:19
  • I just tried with SABS_1314_High.shp and it worked fine – Ian Turton Jan 24 '17 at 14:21
1

You need to reproject your points to be in the SABS projection (US National Atlas Equal Area EPSG:2163):

private boolean isInShape(Point p) throws NoSuchAuthorityCodeException, FactoryException, MismatchedDimensionException, TransformException {


    CoordinateReferenceSystem wgs84 = DefaultGeographicCRS.WGS84;
    CoordinateReferenceSystem layerCRS = features.getSchema().getCoordinateReferenceSystem();
    MathTransform transform= CRS.findMathTransform(wgs84, layerCRS,true);
    Point rp = (Point) JTS.transform(p, transform);
    System.out.println(p+" "+rp);



    Expression propertyName = filterFactory.property(features.getSchema()
            .getGeometryDescriptor().getName());
    Filter filter = filterFactory.contains(propertyName,
            filterFactory.literal(rp));
    SimpleFeatureCollection sub = features.subCollection(filter);
    if (sub.size() > 0) {
        return true;
    }
    return false;
}
  • Thanks, @iant. I understand what you are saying regarding the projections, but I don't know where to get the original points (p) from your example or how I use them once I have them. Should I be creating a new shapefile from the points and using that for searching the points-in-polygons? – Todd Jan 24 '17 at 14:21
  • I'm assuming that Point P in my example is your lat/lon point - then use Point RP to search with – Ian Turton Jan 24 '17 at 14:22
  • Sorry, just clicked. You are suggesting I convert my point for test first, then test it. – Todd Jan 24 '17 at 14:23
  • yes, that's all you need to do – Ian Turton Jan 24 '17 at 14:24
  • that is a simple and concise response that works beautifully, thanks again. Now, I just need to figure out how to determine which polygon the point belongs in and associate that with a school! – Todd Jan 24 '17 at 14:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.