zonebuildings.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)

// The total number of structures in the new maps’ Special Flood Hazard Areas — but not in the old maps’ SPHAs — that were damaged during Sandy.

// basically select * from buildings as b, preliminary_maps as p, existing_county_maps as e where st_intersects(p.the_geom, b.the_geom) and not st_intersects(e.the_geom, b.the_geom) and e.geoid IN ('36085', '36005', '36061', '36081', '36047');
int main() {
  OGRRegisterAll();

  OGRSFDriver *shpDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName("ESRI Shapefile");

  OGRDataSource *out = shpDriver->CreateDataSource("./data/buildings_new_zones.shp", NULL);
  ugh(out, "delete intersections pls.");
  OGRLayer *bldgShpLayer = out->CreateLayer("buildings", NULL, wkbUnknown, NULL);

  OGRDataSource *buildings = OGRSFDriverRegistrar::Open("./data/buildings.shp", FALSE);
  OGRDataSource *pwm = OGRSFDriverRegistrar::Open("./data/pwm.shp", FALSE);
  OGRDataSource *existingMaps = OGRSFDriverRegistrar::Open("./data/eml.shp", FALSE);

  OGRLayer *buildingsLayer = buildings->ExecuteSQL("select * from buildings", NULL, NULL);
  OGRFeatureDefn *bldgFields = buildingsLayer->GetLayerDefn();

  for(int i=0; i < bldgFields->GetFieldCount(); i++){
    OGRFieldDefn *field = bldgFields->GetFieldDefn(i);

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


  OGRLayer *pwmLayer = pwm->ExecuteSQL("select * from pwm", NULL, NULL);
  OGRLayer *existingMapsLayer = existingMaps->ExecuteSQL("select * from eml", NULL, NULL);

  std::cout << buildingsLayer->GetFeatureCount() << std::endl;
  std::cout << pwmLayer->GetFeatureCount() << std::endl;
  std::cout << existingMapsLayer->GetFeatureCount() << std::endl;

  OGRFeature *pwmFeat = pwmLayer->GetFeature(0);
  OGRGeometry *pwmGeom = pwmFeat->GetGeometryRef()->Buffer(0,0);

  OGRFeature *emlFeat = existingMapsLayer->GetFeature(0);
  OGRGeometry *emlGeom = emlFeat->GetGeometryRef();

  //iterate through every building.. to see if it hits the EML and the PWM
  int bldgCt = 0;
  int bldgHitCt = 0;
  int bldg2007Ct = 0;

  OGRFeature *buildingsFeat;
  while((buildingsFeat = buildingsLayer->GetNextFeature()) != NULL) {
    OGRGeometry *buildingsGeom = buildingsFeat->GetGeometryRef();
    if (buildingsGeom->Intersects(emlGeom)) {
      // SKIP buildings that are in existing maps
    } else if (buildingsGeom->Intersects(pwmGeom)) {
      std::cout << "hit " << bldgCt << std::endl;

      // write this feature out to the shp
      OGRGeometry *outGeom;
      outGeom = buildingsGeom;
      OGRFeature *shpFeat = OGRFeature::CreateFeature(bldgFields);


      for(int i=0; i < bldgFields->GetFieldCount(); i++){
        OGRField *bldgField = buildingsFeat->GetRawFieldRef(i);
        OGRFieldDefn *bldgFieldDefn = buildingsFeat->GetFieldDefnRef(i);
        const char *fieldName = bldgFieldDefn->GetNameRef();

        // write out the field
        shpFeat->SetField(i, bldgField);

        const char *field = buildingsFeat->GetFieldAsString(i);
        // std::cout << field << std::endl;
        if (fieldName == std::string("yrb") || fieldName == std::string("yrb_rng") || fieldName == std::string("yra1") || fieldName == std::string("yra1_rng") || fieldName == std::string("yra2") || fieldName == std::string("yra2_rng")) {
          if (atoi(field) > 2006) {
            std::cout << ".......... new bldg " << field << std::endl;
            bldg2007Ct++;
            break;
          }
        }
      }

      shpFeat->SetGeometry(outGeom);
      bldgShpLayer->CreateFeature(shpFeat);
      OGRFeature::DestroyFeature(shpFeat);
      out->SyncToDisk();

      bldgHitCt++;
    }
    bldgCt++;
    OGRFeature::DestroyFeature(buildingsFeat);    
  }
  std::cout << "buildings in preliminary maps, not in existing maps " << bldgHitCt << std::endl;
  std::cout << "buildings in preliminary maps, not in existing maps > 2006 " << bldg2007Ct << std::endl;

}