predictor.cc

#include <ogrsf_frmts.h>
#include <iostream>
#include <fstream>
#include <string>
#include <vector>

#define ugh(it, message) do {\
  if(it == NULL) {\
    std::cout << message << std::endl << std::flush;\
    exit(1);\
  }\
} while(0)


int main() {
  std::string csv = "";
  std::string fips [] = {"36081", "36061", "36047", "36005", "36085", "36059", "36103", "34003", "34017", "34013", "34039", "34023", "34025", "34029", "34005", "34007", "34015", "34033", "34011", "34009", "34001"};
  
  OGRRegisterAll();
  OGRSFDriver *shpDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName("ESRI Shapefile");
  OGRDataSource *sandyZones = OGRSFDriverRegistrar::Open("./data/sandy_zones.shp", FALSE);

  // create the shp
  OGRDataSource *out = shpDriver->CreateDataSource("./data/zone_existing_maps_intersections.shp", NULL);
  ugh(out, "delete intersections pls.");
  OGRLayer *intersects = out->CreateLayer("intersections", NULL, wkbUnknown, NULL);


  for (int i=0; i < 21; i++) {
    std::cout << "DOING: " << fips[i] << std::endl;

    std::string ctyShapefilePath = std::string("./data/unioned/flat_cty_") + fips[i] + std::string(".shp");
    OGRDataSource *existingMaps = OGRSFDriverRegistrar::Open(ctyShapefilePath.c_str(), FALSE);

    std::string sandyGeoidQuery = std::string("SELECT * from sandy_zones where geoid = '") + fips[i] + std::string("'");
    std::string ctyExistingGeoidQuery = std::string("SELECT * from flat_cty_") + fips[i] + std::string("");
    std::cout << ctyExistingGeoidQuery << std::endl;

    OGRLayer *ctySandy = sandyZones->ExecuteSQL(sandyGeoidQuery.c_str(), NULL, NULL);
    OGRLayer *ctyExisting = existingMaps->ExecuteSQL(ctyExistingGeoidQuery.c_str(), NULL, NULL);


    // on first run, create the shapefile layers
    OGRFeatureDefn *fields = ctyExisting->GetLayerDefn();
    if (i == 0) {
      for(int i=0; i < fields->GetFieldCount(); i++){
        OGRFieldDefn *field = fields->GetFieldDefn(i);

        if(intersects->CreateField(field) != OGRERR_NONE) {
          std::cout << "Couldn't make layer" << std::endl; exit(1);
        };
      }
      std::cout << "Copied" << std::endl;
    }


    OGRFeature *ctySandyFeat = ctySandy->GetNextFeature();
    OGRGeometry *ctySandyGeom = ctySandyFeat->GetGeometryRef();
    std::cout << ctySandyGeom->getGeometryType() << std::endl;
    OGRMultiPolygon *ctySandyMGeom = NULL;
    OGRPolygon *ctySandyPGeom = NULL;
    switch(ctySandyGeom->getGeometryType()) {
      case wkbMultiPolygon:
        ctySandyMGeom = (OGRMultiPolygon *) ctySandyGeom;
        break;
      case wkbPolygon:
        ctySandyPGeom = (OGRPolygon *) ctySandyGeom;
        break;
    }

    // area of the county sandy extent
    double ctySandyGeomArea = ctySandyMGeom ? ctySandyMGeom->get_Area() : ctySandyPGeom->get_Area();

    //area of the existing zones' intersection with the sandy extent
    double intersectionArea = 0;
    int featureCt = 0;
    OGRFeature *ctyExistingFeat;
    while((ctyExistingFeat = ctyExisting->GetNextFeature()) != NULL) {
      featureCt++;
      std::cout << featureCt << std::endl;

      OGRMultiPolygon *ctyExistingFeatMGeom = NULL;
      OGRPolygon *ctyExistingFeatPGeom = NULL;
      OGRMultiPolygon *intersectionMGeom = NULL;
      OGRPolygon *intersectionPGeom = NULL;
      OGRGeometry *ctyExistingGeom = ctyExistingFeat->GetGeometryRef();

      //write it to the shp
      OGRGeometry *outGeom;
      outGeom = ctyExistingGeom->Intersection(ctySandyGeom);
      OGRFeature *feat = OGRFeature::CreateFeature(ctyExisting->GetLayerDefn());

      for(int i=0; i < fields->GetFieldCount(); i++){
        OGRField *field = ctyExistingFeat->GetRawFieldRef(i);
        feat->SetField(i, field);
      }

      feat->SetGeometry(outGeom);
      intersects->CreateFeature(feat);
      OGRFeature::DestroyFeature(feat);
      out->SyncToDisk();
     

      switch(ctyExistingGeom->getGeometryType()) {
        case wkbMultiPolygon:
          ctyExistingFeatMGeom = (OGRMultiPolygon *) ctyExistingGeom;
          if (ctySandyMGeom) {
            if (ctyExistingFeatMGeom->Intersects(ctySandyMGeom)) {
              OGRGeometry *intersectionGeom = ctyExistingFeatMGeom->Intersection(ctySandyMGeom);
              if (intersectionGeom->getGeometryType() == wkbMultiPolygon) {
                intersectionMGeom = (OGRMultiPolygon *) intersectionGeom;
                intersectionArea += intersectionMGeom->get_Area();
              } else if (intersectionGeom->getGeometryType() == wkbPolygon) {
                intersectionPGeom = (OGRPolygon *) intersectionGeom;
                intersectionArea += intersectionPGeom->get_Area();
              }
              std::cout << "case1 intersectionArea: " << intersectionArea << std::endl;
            }
          } else {
             if (ctyExistingFeatMGeom->Intersects(ctySandyPGeom)) {
              OGRGeometry *intersectionGeom = ctyExistingFeatMGeom->Intersection(ctySandyPGeom);
              if (intersectionGeom->getGeometryType() == wkbMultiPolygon) {
                intersectionMGeom = (OGRMultiPolygon *) intersectionGeom;
                intersectionArea += intersectionMGeom->get_Area();
              } else if (intersectionGeom->getGeometryType() == wkbPolygon) {
                intersectionPGeom = (OGRPolygon *) intersectionGeom;
                intersectionArea += intersectionPGeom->get_Area();
              }
              std::cout << "case2 intersectionArea: " << intersectionArea << std::endl;
            }
          }
          break;
        case wkbPolygon:
          ctyExistingFeatPGeom = (OGRPolygon *) ctyExistingGeom;
          if (ctySandyMGeom) {
            if (ctyExistingFeatPGeom->Intersects(ctySandyMGeom)) {
              OGRGeometry *intersectionGeom = ctyExistingFeatPGeom->Intersection(ctySandyMGeom);
              if (intersectionGeom->getGeometryType() == wkbMultiPolygon) {
                intersectionMGeom = (OGRMultiPolygon *) intersectionGeom;
                intersectionArea += intersectionMGeom->get_Area();
              } else if (intersectionGeom->getGeometryType() == wkbPolygon) {
                intersectionPGeom = (OGRPolygon *) intersectionGeom;
                intersectionArea += intersectionPGeom->get_Area();
              }
              std::cout << "case3 intersectionArea: " << intersectionArea << std::endl;
            }
          } else {
            if (ctyExistingFeatPGeom->Intersects(ctySandyPGeom)) {
              OGRGeometry *intersectionGeom = ctyExistingFeatPGeom->Intersection(ctySandyPGeom);
              if (intersectionGeom->getGeometryType() == wkbMultiPolygon) {
                intersectionMGeom = (OGRMultiPolygon *) intersectionGeom;
                intersectionArea += intersectionMGeom->get_Area();
              } else if (intersectionGeom->getGeometryType() == wkbPolygon) {
                intersectionPGeom = (OGRPolygon *) intersectionGeom;
                intersectionArea += intersectionPGeom->get_Area();
              }
              std::cout << "case4 intersectionArea: " << intersectionArea << std::endl;
            }
          }
          break;
      }
    }
    std::cout << "===> FIPS: " << fips[i] << std::endl;
    std::cout << "===> EXTENT AREA: " << ctySandyGeomArea << std::endl;
    std::cout << "===> INTERSECTION AREA: " << intersectionArea << std::endl;
    csv += ( fips[i] + std::string(",") + std::to_string(ctySandyGeomArea) + std::string(",") + std::to_string(intersectionArea) + std::string("\n") );
    OGRFeature::DestroyFeature(ctyExistingFeat);
    OGRFeature::DestroyFeature(ctySandyFeat);
  }
  std::ofstream csvFile;
  csvFile.open("./data/predictions.csv");
  csvFile << csv << "\n";
  csvFile.close();
}