用gdal打开一个shp文件并读出里面点的坐标，并计算面积。

#include
#include<stdlib.h>
#include
#include
#include
#include"ogrsf_frmts.h"
using namespace std;
extern double getarea(OGRLinearRing *ring);
int main()
{
GDALAllRegister();
GDALDataset lauds;
lauds = (GDALDataset)GDALOpenEx(“E:/Desktop/national fundamental information/LanduseData/LanduseData/Landuse.shp”, GDAL_OF_VECTOR||GDAL_OF_UPDATE, NULL, NULL, NULL);

if (lauds == NULL)
{
	cout << "open failed\n" << endl;
	exit(1);
}

OGRLayer *laulayer;
laulayer = lauds->GetLayerByName("Landuse");
//laulayer->DeleteField(5);

#if 0
OGRFieldDefn ofield(“area_sum”, OFTReal);
if (laulayer->CreateField(&ofield, TRUE) != OGRERR_NONE)
{
cout << “Creating Name field failed.\n” << endl;
}
#endif

OGRFeature *laufeature;
laulayer->ResetReading();
while ((laufeature=laulayer->GetNextFeature())!=NULL)
{

#if 1
OGRFeatureDefn* laufeaturedf = laulayer->GetLayerDefn();
int ifield;
for (ifield = 0; ifield < laufeaturedf->GetFieldCount(); ifield++)
{
OGRFieldDefn * laufielddf = laufeaturedf->GetFieldDefn(ifield);
if (laufielddf->GetType() == OFTInteger)
cout << laufielddf->GetNameRef() <<":"<< laufeature->GetFieldAsInteger(ifield) << endl;
else if (laufielddf->GetType() == OFTReal)
cout << laufielddf->GetNameRef() <<":"<< laufeature->GetFieldAsDouble(ifield) << endl;
else if (laufielddf->GetType() == OFTString)
cout << laufielddf->GetNameRef() << “:”<GetFieldAsString(ifield) << endl;
else
cout << “else:” << laufeature->GetFieldAsString(ifield);
}
#endif

	double area=0;
	OGRGeometry *laugeometry;
	laugeometry = laufeature->GetGeometryRef();
	if (laugeometry != NULL)
	{
		OGRwkbGeometryType laugeotype = laugeometry->getGeometryType();
		if (laugeotype == wkbPoint)
		{
			OGRPoint *laupoint = (OGRPoint *)laugeometry;
			cout << laupoint->getX() << laupoint->getY() << endl;
		}

		else if (laugeotype == wkbPolygon)
		{
			cout << "polygon" << endl;
			OGRPolygon *laupolygon = (OGRPolygon *)laugeometry->clone();
			OGRLinearRing  *laulr = laupolygon->getExteriorRing();
			area += getarea(laulr);


			if (laupolygon->getNumInteriorRings()) {
				for (int i = 0; i < laupolygon->getNumInteriorRings(); i++)
				{
					OGRLinearRing *laulr = (OGRLinearRing*)laupolygon->getInteriorRing(i);
					area += getarea(laulr);
				}
			}


		}

		else if(laugeotype == wkbMultiPolygon)
		{
			cout << "multipolygon" << endl;
			OGRMultiPolygon *laumltupolygon = (OGRMultiPolygon *)laugeometry;
			OGRPolygon *laupolygon1 = NULL;
			for (int i = 0; i < laumltupolygon->getNumGeometries(); i++)
			{
				laupolygon1 = (OGRPolygon*)laumltupolygon->getGeometryRef(i);
				OGRLinearRing *laulr = (OGRLinearRing*)laupolygon1->getExteriorRing();
				area += getarea(laulr);
				if (laupolygon1->getNumInteriorRings()) {
					for(int i=0;i<laupolygon1->getNumInteriorRings();i++)
					{
						OGRLinearRing *laulr = (OGRLinearRing*)laupolygon1->getInteriorRing(i);
						area += getarea(laulr);
					}
				}
			}
		}
		else cout << "no geometry" << endl;



	}

	
	cout << "areaaaa:\t"<<area<<"\n"<<"end" << "\n"<<endl;

#if 0
laufeature->SetField(“area_sum”, area); // Set fields
laulayer->SetFeature(laufeature);
#endif

	OGRFeature::DestroyFeature(laufeature);
}
cout << "all end" << endl;
GDALClose(lauds);
system("pause");

}

gearea。cpp
#include
#include
#include"ogrsf_frmts.h"
using namespace std;
double getarea(OGRLinearRing ring)
{
double X1, X2, Y1, Y2;
double pi = 3.14159265358979, rou = 206264.8062471, lr = 6378137.0, sr = 6356752.314, one,/=6.6943800229008E-03,*/two = 6.7394967754790E-03, ji = 6399593.6258640232;
double afa = 1 / 298.257222101;//afa-椭球扁率；；pi-π；rou=p；lr=长半轴；sr=短半轴；one=椭球第一偏心率；two=椭球第二偏心率；ji=极点子午圈曲率半径；
double A, B, C, D, E, lad_, lad, log, log_, s = 0;
one = (pow(lr, 2) - pow(sr, 2)) / (pow(lr, 2));
A = 1 + (3.0 / 6.0)*one + (30.0 / 80.0)*pow(one, 2.0) + (35.0 / 112.0)*pow(one, 3) + (630.0 / 2304.0)*pow(one, 4.0);
B = (1 / 6.0)*one + (15 / 80.0)*pow(one, 2) + (21.0 / 112)*pow(one, 3) + (420.0 / 2304)*pow(one, 4);
C = (3.0 / 80)*pow(one, 2) + (7.0 / 112)*pow(one, 3) + (180.0 / 2304)*pow(one, 4);
D = (1.0 / 112)*pow(one, 3) + (45.0 / 2304)*pow(one, 4);
E = (5.0 / 2304)*pow(one, 4);
int n, i = ring->getNumPoints();

for (n = 0; n < i - 1; n++)
{

	X1 = ring->getX(n + 1)*pi / 180.;
	X2 = ring->getX(n)*pi / 180.;
	Y1 = ring->getY(n + 1)*pi / 180.;
	Y2 = ring->getY(n)*pi / 180.;
	lad = (X2 + X1) / 2.;
	lad_ = 2. * sr*sr*lad;
	log_ = (Y2 - Y1) / 2.;
	log = (Y2 + Y1) / 2.;
	s = s + lad_ * (A*sin(log_)*cos(log) - B * sin(3 * log_)*cos(3 * log) + C * sin(5 * log_)*cos(5 * log) - D * sin(7 * log_)*cos(7 * log) + E * sin(9 * log_)*cos(9 * log));

}
if (n = i - 1)
{
	X1 = ring->getX(0)*pi / 180.;
	X2 = ring->getX(i - 1)*pi / 180.;
	Y1 = ring->getY(0)*pi / 180.;
	Y2 = ring->getY(i - 1)*pi / 180.;
	lad = (X2 + X1) / 2;
	lad_ = 2 * sr*sr*lad;
	log_ = (Y2 - Y1) / 2;
	log = (Y2 + Y1) / 2;
	s += lad_ * (A*sin(log_)*cos(log) - B * sin(3 * log_)*cos(3 * log) + C * sin(5 * log_)*cos(5 * log) - D * sin(7 * log_)*cos(7 * log) + E * sin(9 * log_)*cos(9 * log));

}

return s;

}