GEOTIFF图片切割

#include "stdafx.h" #pragma comment(lib, "gdal_i.lib") #include "include/gdal/gdal.h" #include "include/gdal/gdal_priv.h" #include "iostream" #include "fstream" #include <io.h> #include <string.h> #include <stdio.h> #include <errno.h> #include "include/gdal/ogr_spatialref.h" #include "include/gdal/cpl_minixml.h" #include "include/gdal/gdal_alg.h" #include "include/gdal/gdalwarper.h" #include <afx.h> #include "math.h" #define PI 3.1415926535 #define originShift 20037508.342789244 #define initialResolution 156543.03392804097 #define tileSize 256 using namespace std; typedef unsigned char BYTE; double cost_time; clock_t start,end; class Mercator { public: double mx,my;//墨卡托 double lat,lon;//经纬度 double tx,ty;//瓦片 double px,py;//像素 double zoom;//级数 double px6,py6; double minx,miny;//最小点 double maxx,maxy;//最大点 double minlat,minlon;//第9个函数 double maxlat,maxlon; public: void LatLonToMeters(double lat,double lon); void MetersToLatLon(double mx,double my ); void PixelsToMeters(double px,double py,double zoom); void MetersToPixels(double mx,double my,double zoom); void PixelsToTile(double px,double py); void PixelsToRaster(double px,double py,double zoom); void MetersToTile(double mx,double my,double zoom); void TileBounds(double tx,double ty,double zoom); void TileLatLonBounds(double tx,double ty,double zoom ); double Resolution(double zoom ); int ZoomForPixelSize(double pixelSize ); protected: private: }; void Mercator::LatLonToMeters(double lat,double lon) { double mxx= lon *originShift/180.0; double myy= log(tan((90 +lat) * PI/ 360.0 ))/(PI/180.0); myy = myy *originShift/180.0; mx=mxx; my=myy; } void Mercator::MetersToLatLon(double mx,double my) { double lonn = (mx / originShift) * 180.0; double latt = (my / originShift) * 180.0; latt = 180 /PI * (2 * atan(exp( latt * PI / 180.0)) -PI / 2.0); lon=lonn; lat=latt; } void Mercator::PixelsToMeters(double px,double py,double zoom) { double res=initialResolution /(pow(2,zoom)); double mxx = px * res - originShift; double myy = py * res - originShift; mx=mxx; my=myy; } void Mercator::MetersToPixels(double mx,double my,double zoom) { double res=initialResolution /(pow(2,zoom)); double pxx = (mx + originShift) / res; double pyy = (my + originShift) / res; px=pxx; py=pyy; } void Mercator::PixelsToTile(double px,double py) { double txx = int( ceil( px / float(tileSize) ) - 1 ); double tyy = int( ceil( py / float(tileSize) ) - 1 ); tx=txx; ty=tyy; } void Mercator::PixelsToRaster(double px,double py,double zoom) { double mapSize = tileSize*(pow(2,zoom)); px6=px; py6=mapSize - py; } void Mercator::MetersToTile(double mx,double my,double zoom) { double res=initialResolution /(pow(2,zoom)); double pxx=(mx+originShift) / res; double pyy=(my+originShift) / res; px=pxx; py=pyy; double txx = int( ceil( px / float(tileSize) ) - 1 ); double tyy = int( ceil( py / float(tileSize) ) - 1 ); tx=txx; ty=tyy; } void Mercator::TileBounds(double tx,double ty,double zoom) { //minx, miny = PixelsToMeters( tx*tileSize, ty*tileSize, zoom ) double res=initialResolution /(pow(2,zoom)); double mxx = tx*tileSize * res - originShift; double myy = ty*tileSize * res - originShift; minx=mxx; miny=myy; //maxx, maxy = PixelsToMeters( (tx+1)*tileSize, (ty+1)*tileSize, zoom ) //double res=initialResolution /(pow(2,zoom)); mxx = (tx+1)*tileSize* res - originShift; myy = (ty+1)*tileSize* res - originShift; maxx=mxx; maxy=myy; } void Mercator::TileLatLonBounds(double tx,double ty,double zoom ) { //bounds = TileBounds( tx, ty, zoom) double res=initialResolution /(pow(2,zoom)); double mxx = tx*tileSize * res - originShift; double myy = ty*tileSize * res - originShift; minx=mxx; miny=myy; mxx = (tx+1)*tileSize* res - originShift; myy = (ty+1)*tileSize* res - originShift; maxx=mxx; maxy=myy; //minLat, minLon = MetersToLatLon(bounds[0], bounds[1]) double lonn = (minx / originShift) * 180.0; double latt = (miny / originShift) * 180.0; latt = 180 /PI * (2 * atan(exp( lat * PI / 180.0)) -PI / 2.0); minlon=lonn; minlat=latt; //maxLat, maxLon = MetersToLatLon(bounds[2], bounds[3]) lonn = (maxx / originShift) * 180.0; latt = (maxy / originShift) * 180.0; latt = 180 /PI * (2 * atan(exp( lat * PI / 180.0)) -PI / 2.0); maxlon=lonn; maxlat=latt; } double Mercator::Resolution(double zoom ) { return initialResolution /(pow(2,zoom)); } int Mercator::ZoomForPixelSize(double pixelSize ) { for (int i=0;i<18;i++) { if (pixelSize>Resolution(i)) { if(i!=0) return i-1; } //else // /*return 0;*/ } } struct Zdzx { int rx, ry, rxsize, rysize; int wx, wy, wxsize, wysize; }; void geo_query( double minx, double maxy,double miny,double maxx ,double omaxy,int xsize,int ysize,double xx,double xx1,double xx5,Zdzx &cx) { int rx, ry, rxsize, rysize; int wx, wy, wxsize, wysize; int rxshift=0,ryshift=0; wxsize = 256; wysize = 256; rx= int((minx - xx)/xx1 + 0.001); ry= int((maxy - omaxy)/xx5 + 0.001); rxsize= (int)( (maxx - minx) / xx1+ 0.5); rysize= (int)( (miny - maxy) / xx5+ 0.5); wx = 0; if (rx<0) { rxshift = abs(rx); wx = (int)( wxsize * rxshift / rxsize+0.001) ; wxsize = wxsize - wx; rxsize = rxsize -( int)( rxsize * rxshift / rxsize+0.5) ; rx=0; }; if (rx+rxsize>xsize) { wxsize = (int)( wxsize * (xsize - rx) / rxsize); rxsize = xsize - rx; }; wy=0; if (ry<0) { ryshift = abs(ry); wy = (int)( wysize * ryshift/ rysize ); wysize = wysize - wy; rysize = rysize -( int)( rysize * ryshift/ rysize ); ry = 0; }; if (ry+rysize > ysize) { wysize = (int)( wysize * (ysize - ry) / rysize); rysize = ysize - ry; }; cx.rx=rx; cx.ry=ry; cx.rxsize=rxsize; cx.rysize=rysize; cx.wx=wx; cx.wxsize=wxsize; cx.wy=wy; cx.wysize=wysize; }; CString CreateFolder(CString path1,CString path2,CString path3,CString filename,CString fileExtension) { CString path=path1; path+="//"+path2; CreateDirectory(path,NULL); path+="//"+path3; CreateDirectory(path,NULL); path+="//"+filename+fileExtension; return path; }; int main() { CPLErr errocode;//错误标记 start=clock(); //开始的时候记录时间 GDALAllRegister(); GDALDataset *poDataset; poDataset = (GDALDataset *) GDALOpen( "d:/test/jl1.tif", GA_ReadOnly ); const char * oTargetSRS=NULL; const char * oSourceSRS=GDALGetProjectionRef(poDataset); CString ss; ss="PROJCS[\"Google Maps Global Mercator\",\n\ GEOGCS[\"WGS 84\",\n\ DATUM[\"WGS_1984\",\n\ SPHEROID[\"WGS 84\",6378137,298.2572235630016,\n\ AUTHORITY[\"EPSG\",\"7030\"]],\n\ AUTHORITY[\"EPSG\",\"6326\"]],\n\ PRIMEM[\"Greenwich\",0],\n\ UNIT[\"degree\",0.0174532925199433],\n\ AUTHORITY[\"EPSG\",\"4326\"]],\n\ PROJECTION[\"Mercator_1SP\"],\n\ PARAMETER[\"central_meridian\",0],\n\ PARAMETER[\"scale_factor\",1],\n\ PARAMETER[\"false_easting\",0],\n\ PARAMETER[\"false_northing\",0],\n\ UNIT[\"metre\",1,\n\ AUTHORITY[\"EPSG\",\"9001\"]]]"; oTargetSRS=ss; GDALDataset *poDataset5= (GDALDataset *)GDALAutoCreateWarpedVRT ( poDataset, oSourceSRS, oTargetSRS, GRA_NearestNeighbour, 0.0, NULL); double adfGeoTransform[6],zuichude[6]; poDataset5->GetGeoTransform( adfGeoTransform ) ; poDataset->GetGeoTransform( zuichude ) ; /////////////////////////////////////////////////////////////////////////// double ominx,ominy,omaxx,omaxy; int xsize,ysize,tilesize=256; xsize=poDataset5->GetRasterXSize(); ysize=poDataset5->GetRasterYSize(); ominx = adfGeoTransform[0]; omaxx = adfGeoTransform[0]+xsize*adfGeoTransform[1]; omaxy= log(tan((90 +zuichude[3]) * PI/ 360.0 ))/(PI/180.0)*o