[Gdal-dev] re: How To Specify GCP for gdalwarp
Matt Lynch
matt at terraEngine.com
Wed Aug 13 19:24:36 EDT 2003
Hi Bill,
Since the listserv seems so quiet, I will take a stab at it (disclaimer
- I'm still figuring this out). My read of the comments on the gdalwarp
utility make me think that you need to have the GCP already assigned to
the image. The GDAL API does have routines to add GCPs to an image.
I took an initial stab at doing this for you, but ran into trouble with
the JPEG driver. Perhaps you could convert to a tif before running?
The run is here:
C:\matt\gdal-cvs-030403\apps>addGCP a.jpg -120.6299 56.3091 1 1
-119.7047 56.145 9 500 1 -120.0005 55.6328 500 500 -120.9143 55.7944 1.0
500 -120.3118 55.9714 25 0 250 setting the following ground control
points for a.jpg:
gcp0 (-120.629900,56.309100) tied to (1.000000,1.000000).
gcp1 (-119.704700,56.145900) tied to (500.000000,1.000000).
gcp2 (-120.000500,55.632800) tied to (500.000000,500.000000).
gcp3 (-120.914300,55.794400) tied to (1.000000,500.000000).
gcp4 (-120.311800,55.971400) tied to (250.000000,250.000000).
ERROR 6: The JPEG driver does not support update access to existing
datasets.
C:\matt\gdal-cvs-030403\apps>addGCP a.tif -120.6299 56.3091 1 1
-119.7047 56.145
9 500 1 -120.0005 55.6328 500 500 -120.9143 55.7944 1.0 500 -120.3118
55.9714 25
0 250
setting the following ground control points for a.tif:
gcp0 (-120.629900,56.309100) tied to (1.000000,1.000000).
gcp1 (-119.704700,56.145900) tied to (500.000000,1.000000).
gcp2 (-120.000500,55.632800) tied to (500.000000,500.000000).
gcp3 (-120.914300,55.794400) tied to (1.000000,500.000000).
gcp4 (-120.311800,55.971400) tied to (250.000000,250.000000).
ERROR 1: Unable to compute a transformation between pixel/line
and georeferenced coordinates for a.tif.
There is no affine transformation and no GCPs.
ERROR 5: panSrcBands[0] = 3080592 ... out of range for dataset.
Hopefully Frank or others can correct me if I am wrong on gdalwarp.exe
adding GCPs, and maybe clean the code some.
The code is also attached to make it easier for you to copy.
Best of luck,
Matt
#include <stdio.h>
#include <stdlib.h>
#include "gdalwarper.h"
#include "cpl_string.h"
int main(int argc, char *argv[])
{
int i,j;
char *pszSourceFileName = NULL;
GDAL_GCP *gcps = NULL;
GDALDatasetH hSrcDS,hDstDS;
GDALWarpOptions *psWarpOptions;
GDALWarpOperation oOperation;
if( (argc > 2)&&(argc%4==2) )
{
// request memory to hold the GCPs
gcps = (GDAL_GCP *) malloc(
((argc-2)/4)*sizeof(GDAL_GCP) );
pszSourceFileName = argv[1];
for( i=2,j=0 ; i<argc ; i+=4,j++ )
{
gcps[j].dfGCPPixel = atof(argv[i]);
gcps[j].dfGCPLine = atof(argv[i+1]);
gcps[j].dfGCPX = atof(argv[i+2]);
gcps[j].dfGCPY = atof(argv[i+3]);
gcps[j].pszId = (char *) malloc( 20*
sizeof(char) );
sprintf(gcps[j].pszId,"gcp%i",j);
}
}
fprintf(stdout,"setting the following ground control points for
%s:\n",argv[1] );
for(i=0 ; i<j ; i++)
{
fprintf(stdout,"%s\t(%lf,%lf) tied to
(%lf,%lf).\n",gcps[i].pszId,
gcps[i].dfGCPPixel,
gcps[i].dfGCPLine,
gcps[i].dfGCPX,
gcps[i].dfGCPY);
}
GDALAllRegister();
hSrcDS = GDALOpen(pszSourceFileName,GA_Update);
GDALSetGCPs ( hSrcDS, j, gcps, "+proj=latlong +datum=WGS84" ) ;
hDstDS = GDALOpen("AB.tif",GA_Update);
psWarpOptions = GDALCreateWarpOptions();
psWarpOptions->hSrcDS = hSrcDS;
psWarpOptions->hDstDS = hDstDS;
psWarpOptions->nBandCount = 1;
psWarpOptions->panSrcBands = (int *) CPLMalloc(
sizeof(int)*psWarpOptions->nBandCount );
psWarpOptions->panDstBands = (int *) CPLMalloc(
sizeof(int)*psWarpOptions->nBandCount );
psWarpOptions->panDstBands[0] = 1;
psWarpOptions->pfnProgress = GDALTermProgress;
psWarpOptions->pTransformerArg =
GDALCreateGenImgProjTransformer(
hSrcDS,GDALGetProjectionRef(hSrcDS),hDstDS,GDALGetProjectionRef(hDstDS),
FALSE,0.0,1);
psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
oOperation.Initialize(psWarpOptions);
oOperation.ChunkAndWarpImage(0,0,GDALGetRasterXSize(hDstDS),GDALGetRaste
rYSize(hDstDS));
GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg);
GDALDestroyWarpOptions( psWarpOptions );
GDALClose( hDstDS );
GDALClose( hSrcDS );
return 0;
}
-----Original Message-----
From: gdal-dev-admin at remotesensing.org
[mailto:gdal-dev-admin at remotesensing.org] On Behalf Of Bill Bob
Sent: Wednesday, August 13, 2003 10:08 AM
To: gdal-dev at remotesensing.org
Subject: [Gdal-dev] re: How To Specify GCP for gdalwarp
Hi,
I am still trying to figure out the syntax to do this in gdalwarp. In
my
previous question I forgot to mention I was using the dos version of the
utility. What I am not clear on is where and how to specify the ground
control points. Do I do it with the -s_srs parameter or perhaps as a
-co
parameter?
I have searched the gdal website thoroughly, but can't seem to get the
command line right.
Any help is most appreciated!
BB
_________________________________________________________________
MSN 8 with e-mail virus protection service: 2 months FREE*
http://join.msn.com/?page=features/virus
_______________________________________________
Gdal-dev mailing list
Gdal-dev at remotesensing.org
http://remotesensing.org/mailman/listinfo/gdal-dev
---
Incoming mail is certified Virus Free.
Checked by AVG anti-virus system (http://www.grisoft.com).
Version: 6.0.507 / Virus Database: 304 - Release Date: 8/4/2003
---
Outgoing mail is certified Virus Free.
Checked by AVG anti-virus system (http://www.grisoft.com).
Version: 6.0.507 / Virus Database: 304 - Release Date: 8/4/2003
-------------- next part --------------
#include <stdio.h>
#include <stdlib.h>
#include "gdalwarper.h"
#include "cpl_string.h"
int main(int argc, char *argv[])
{
int i,j;
char *pszSourceFileName = NULL;
GDAL_GCP *gcps = NULL;
GDALDatasetH hSrcDS,hDstDS;
GDALWarpOptions *psWarpOptions;
GDALWarpOperation oOperation;
if( (argc > 2)&&(argc%4==2) )
{
// request memory to hold the GCPs
gcps = (GDAL_GCP *) malloc( ((argc-2)/4)*sizeof(GDAL_GCP) );
pszSourceFileName = argv[1];
for( i=2,j=0 ; i<argc ; i+=4,j++ )
{
gcps[j].dfGCPPixel = atof(argv[i]);
gcps[j].dfGCPLine = atof(argv[i+1]);
gcps[j].dfGCPX = atof(argv[i+2]);
gcps[j].dfGCPY = atof(argv[i+3]);
gcps[j].pszId = (char *) malloc( 20* sizeof(char) );
sprintf(gcps[j].pszId,"gcp%i",j);
}
}
fprintf(stdout,"setting the following ground control points for %s:\n",argv[1] );
for(i=0 ; i<j ; i++)
{
fprintf(stdout,"%s\t(%lf,%lf) tied to (%lf,%lf).\n",gcps[i].pszId,
gcps[i].dfGCPPixel,
gcps[i].dfGCPLine,
gcps[i].dfGCPX,
gcps[i].dfGCPY);
}
GDALAllRegister();
/*
CPLErr GDALSetGCPs ( GDALDatasetH hDS, int nGCPCount, const GDAL_GCP * pasGCPList, const char * pszGCPProjection ) ;
*/
hSrcDS = GDALOpen(pszSourceFileName,GA_Update);
GDALSetGCPs ( hSrcDS, j, gcps, "+proj=latlong +datum=WGS84" ) ;
hDstDS = GDALOpen("AB.tif",GA_Update);
psWarpOptions = GDALCreateWarpOptions();
psWarpOptions->hSrcDS = hSrcDS;
psWarpOptions->hDstDS = hDstDS;
psWarpOptions->nBandCount = 1;
psWarpOptions->panSrcBands = (int *) CPLMalloc( sizeof(int)*psWarpOptions->nBandCount );
psWarpOptions->panDstBands = (int *) CPLMalloc( sizeof(int)*psWarpOptions->nBandCount );
psWarpOptions->panDstBands[0] = 1;
psWarpOptions->pfnProgress = GDALTermProgress;
psWarpOptions->pTransformerArg = GDALCreateGenImgProjTransformer( hSrcDS,
GDALGetProjectionRef(hSrcDS),
hDstDS,
GDALGetProjectionRef(hDstDS),
FALSE,
0.0,1);
psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
oOperation.Initialize(psWarpOptions);
oOperation.ChunkAndWarpImage(0,0,GDALGetRasterXSize(hDstDS),GDALGetRasterYSize(hDstDS));
GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg);
GDALDestroyWarpOptions( psWarpOptions );
GDALClose( hDstDS );
GDALClose( hSrcDS );
return 0;
}
More information about the Gdal-dev
mailing list