Steve, The -te option mentions the extents of the output file, not the input file. As of now there is no option in gdalwrap to specify the input file extents.
Regards, -- Chaitanya kumar CH. On Thu, Mar 19, 2009 at 10:44 PM, Stephen Woodbridge <wood...@swoodbridge.com> wrote: > Hi all, > > This is my annual wake-up got to do something with rasters and forgot the > finer details of the last time I did a problem similar to it. > > I have a bunch of LandSat mrsid imagery, they are in multiple UTM zone > projections and I want to convert them to geographic projection, WGS84, and > into tiff files. The files are too large to convert into a single tif, so I > plan to quarter each file into 4 pieces. > > So I wrote a script the reads the files and collects info about them, then > generates 4 gdalwarp commands to reproject each quarter of the file (perl > included below). A typical command looks like: > > gdalwarp -srcnodata 0 -dstnodata 0 -s_srs +init=epsg:32642 -t_srs EPSG:4326 > -te 68.9999676729652 24.9707453133203 72.7715338115736 27.4799559407707 -rb > -wm 250 --config GDAL_ONE_BIG_READ ON -co "TILED=YES" > ./GeoCover2000/n-42-25.sid ./af/n-42-25_1-0.tif > > > The reprojection seems to have worked fine, but the bounds seem to be off > because of the black (nodata) areas. > > You can see my poor results here: > http://imaptools.com/downloads/raster.gif > > Is there a better way to do this? I seem to remember a tool to generate > tiles. > Does this need to be done in two steps? How? > Did I mess up my logic below? I've been over it multiple times and do not > see any errors. > > Any help would be appreciated. > > Thanks, > -Steve > > $ gdalinfo -nomd ./GeoCover2000/n-42-25.sid > Driver: MrSID/Multi-resolution Seamless Image Database (MrSID) > Files: ./GeoCover2000/n-42-25.sid > Size is 51078, 39090 > Coordinate System is `' > Origin = (136066.125000000000000,3323577.375000000000000) > Pixel Size = (14.250000000000000,-14.250000000000000) > Corner Coordinates: > Upper Left ( 136066.125, 3323577.375) > Lower Left ( 136066.125, 2766544.875) > Upper Right ( 863927.625, 3323577.375) > Lower Right ( 863927.625, 2766544.875) > Center ( 499996.875, 3045061.125) > Band 1 Block=1024x128 Type=Byte, ColorInterp=Red > Minimum=0.000, Maximum=226.000, Mean=109.610, StdDev=39.202 > Overviews: 25539x19545, 12770x9773, 6385x4887, 3193x2444, 1597x1222, > 799x611, 400x306, 200x153 > Band 2 Block=1024x128 Type=Byte, ColorInterp=Green > Minimum=0.000, Maximum=229.000, Mean=125.833, StdDev=31.828 > Overviews: 25539x19545, 12770x9773, 6385x4887, 3193x2444, 1597x1222, > 799x611, 400x306, 200x153 > Band 3 Block=1024x128 Type=Byte, ColorInterp=Blue > Minimum=0.000, Maximum=227.000, Mean=107.092, StdDev=37.230 > Overviews: 25539x19545, 12770x9773, 6385x4887, 3193x2444, 1597x1222, > 799x611, 400x306, 200x153 > > > $ cat ./GeoCover2000/n-42-25.prj > Projection UTM > Datum WGS84 > Zunits NO > Units METERS > Zone 42 > Xshift 0.000000 > Yshift 0.000000 > Parameters > > > > Perl snippet: > > foreach my $sid (@files) { > > next unless $wanted && $sid =~ /$wanted/i; > > if (! -r $sid) { > warn "SKIPPING: $sid is not readable!\n"; > next; > } > > my $gdalinfo = `gdalinfo -nomd $sid`; > > # parse out the strings that we need > > $gdalinfo =~ /Upper Left \((\s*[^,]+),\s*([^\)]+)/s || warn "Upper > Left parse error\n"; > my $ul = [$1, $2]; > $gdalinfo =~ /Lower Left \((\s*[^,]+),\s*([^\)]+)/s || warn "Lower > Left parse error\n"; > my $ll = [$1, $2]; > $gdalinfo =~ /Upper Right \((\s*[^,]+),\s*([^\)]+)/s || warn "Upper > Right error\n"; > my $ur = [$1, $2]; > $gdalinfo =~ /Lower Right \((\s*[^,]+),\s*([^\)]+)/s || warn "Lower > Right error\n"; > my $lr = [$1, $2]; > > my $outf = $sid; > $outf =~ s/\.sid//; > $outf =~ m/\/([^\/]+)$/; > $outf = $1; > > my $prj = $sid; > $prj =~ s/\.sid/.prj/; > if (! -r $prj) { > warn "SKIPPING: $prj is not readable!\n"; > next; > } > my $projinfo = `grep Zone $prj`; > if ($projinfo =~ /^Zone\s*([-+]?\d+)/) { > my $znum = $1; > if ($znum < 0) { > $epsg = 32700 - $znum; > } else { > $epsg = 32600 + $znum; > } > } > else { > warn "SKIPPING: could not find zone or epsg in $prj\n"; > next; > } > > my $proj = sprintf("+init=epsg:%d", $epsg); > print "Using projection: $proj for $projinfo" if $DEBUG; > > my $pj = Geo::Proj4->new($proj); > if (! $pj) { > warn "SKIPPING: Geo::Proj4 error: " . Geo::Proj4->error . "\n"; > next; > } > > $ul = [reverse $pj->inverse(@$ul)]; > $ll = [reverse $pj->inverse(@$ll)]; > $ur = [reverse $pj->inverse(@$ur)]; > $lr = [reverse $pj->inverse(@$lr)]; > > # this is probably a parallelgram in latlong > if ($DEBUG) { > printf("ul = %11.6f, %11.6f\n", @$ul); > printf("ll = %11.6f, %11.6f\n", @$ll); > printf("ur = %11.6f, %11.6f\n", @$ur); > printf("lr = %11.6f, %11.6f\n", @$lr); > } > > # Get the bbox of the parallelgram > my $xmin = $ul->[0]<$ll->[0]?$ul->[0]:$ll->[0]; > my $xmax = $ur->[0]>$lr->[0]?$ur->[0]:$lr->[0]; > my $ymin = $ll->[1]<$lr->[1]?$ll->[1]:$lr->[1]; > my $ymax = $ul->[1]>$ur->[1]?$ul->[1]:$ur->[1]; > > my $tiled = '-co "TILED=YES"'; > $tiled = '' if -e $outf; > > my $nparts = 2; > my $dx = ($xmax - $xmin) / $nparts; > my $dy = ($ymax - $ymin) / $nparts; > my ($x0, $x1, $y0, $y1); > for (my $i=0, $x0=$xmin; $i<$nparts; $i++, $x0 += $dx) { > $x1 = $x0 + $dx; > for (my $j=0, $y0=$ymin; $j<$nparts; $j++, $y0 += $dy) { > $y1 = $y0 + $dy; > my $tif = sprintf("%s/%s/%s_%d-%d.tif", $work, $outd, $outf, > $i, $j); > my $cmd = sprintf("gdalwarp -srcnodata 0 -dstnodata 0 -s_srs > +init=epsg:%d -t_srs EPSG:4326 -te $x0 $y0 $x1 $y1 -rb -wm 250 --config > GDAL_ONE_BIG_READ ON $tiled $sid $tif\n", $epsg); > > print "$cmd\n"; > my $rc = system($cmd); > print "SYSTEM returned: " . sprintf("%lx\n", ($rc & 0xffff)) > if $DEBUG > } > } > } > > _______________________________________________ > gdal-dev mailing list > gdal-dev@lists.osgeo.org > http://lists.osgeo.org/mailman/listinfo/gdal-dev > _______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
