Skip to content

Instantly share code, notes, and snippets.

@ashaw
Last active Dec 27, 2015
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)
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();
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment