Skip to content

Instantly share code, notes, and snippets.

@ashaw

ashaw/zonebuildings.cc Secret

Created Nov 22, 2013
Embed
What would you like to do?
#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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment