[gdal-dev] Handling data raster data in decimal seconds

Even Rouault even.rouault at mines-paris.org
Fri Nov 22 10:38:14 PST 2013


Le vendredi 22 novembre 2013 13:12:32, Moses.Gone at t-systems.com a écrit :
> Hello GDAL Experts,
> We have a grid dataset that is in WGS84 decimal-second coordinates. The
> first problem is that when we run gdalinfo against this data, the
> coordinates systems comes out wrong i.e in decimal degrees rather than in
> decimal seconds:
> 
> bash-2.05$ gdalinfo -noct
> /users/ds57201a/user3/rbda/de/uka30qs/dlu/grd/aa_ds_1 Driver: AIG/Arc/Info
> Binary Grid
> Files: /users/ds57201a/user3/rbda/de/uka30qs/dlu/grd/aa_ds_1
>        /users/ds57201a/user3/rbda/de/uka30qs/dlu/grd/aa_ds_1.aux.xml
>        /users/ds57201a/user3/rbda/de/uka30qs/dlu/grd/aa_ds_1/hdr.adf
>        /users/ds57201a/user3/rbda/de/uka30qs/dlu/grd/aa_ds_1/metadata.xml
>        /users/ds57201a/user3/rbda/de/uka30qs/dlu/grd/aa_ds_1/prj.adf
>        /users/ds57201a/user3/rbda/de/uka30qs/dlu/grd/aa_ds_1/w001001.adf
>        /users/ds57201a/user3/rbda/de/uka30qs/dlu/grd/aa_ds_1/w001001x.adf
>        /users/ds57201a/user3/rbda/de/uka30qs/dlu/grd/aa_ds_1/dblbnd.adf
>        /users/ds57201a/user3/rbda/de/uka30qs/dlu/grd/aa_ds_1/sta.adf
>        /users/ds57201a/user3/rbda/de/uka30qs/dlu/grd/aa_ds_1/vat.adf
> Size is 1324, 936
> Coordinate System is:
> GEOGCS["WGS 84",
>     DATUM["WGS_1984",
>         SPHEROID["WGS 84",6378137,298.257223563,
>             AUTHORITY["EPSG","7030"]],
>         TOWGS84[0,0,0,0,0,0,0],
>         AUTHORITY["EPSG","6326"]],
>     PRIMEM["Greenwich",0,
>         AUTHORITY["EPSG","8901"]],
>     UNIT["degree",0.0174532925199433,
>         AUTHORITY["EPSG","9108"]],
>     AUTHORITY["EPSG","4326"]]
> Origin = (5.984444444444445,50.945000000000000)
> Pixel Size = (0.000277777777778,-0.000277777777778)
> Corner Coordinates:
> Upper Left  (   5.9844444,  50.9450000) (  5d59' 4.00"E, 50d56'42.00"N)
> Lower Left  (   5.9844444,  50.6850000) (  5d59' 4.00"E, 50d41' 6.00"N)
> Upper Right (   6.3522222,  50.9450000) (  6d21' 8.00"E, 50d56'42.00"N)
> Lower Right (   6.3522222,  50.6850000) (  6d21' 8.00"E, 50d41' 6.00"N)
> Center      (   6.1683333,  50.8150000) (  6d10' 6.00"E, 50d48'54.00"N)
> 
> Source SRS should look like this:
> GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS
> 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326
> "]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["arc-second",0.0000
> 04848136811095,AUTHORITY["EPSG","9104"]],AUTHORITY["EPSG","4326"]]
> 
> We know for sure that the data is in DS and not DD

I can see in the AIGRID driver that there's a special case to deal with "Units 
DS" when reading in prj.adf.
I'm not sure when your source SRS comes from ?
It would be interesting that you create a ticket with the dataset attached and 
the expect corner coordinates.


> 
> Now the real problem comes when we want to warp this dataset to a different
> SRS. Since GDAL already reads a wrong SRS on the input data, all
> transformations end up wrong. We tried to override this by forcing a
> source coordinate systems on the data using gdal_translate, but to no
> avail:
> 
> gdal_translate -of vrt -a_srs
> /users/ds57201a/user1/uka/ukatiler/config/projections/ds_4326.prj
> aachen.vrt test.vrt
> 
> ds_4326.prj wkt:
> GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS
> 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326
> "]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["arc-second",0.0000
> 04848136811095,AUTHORITY["EPSG","9104"]],AUTHORITY["EPSG","4326"]]
> 
> gdalwarp --debug on -wm 500 -wo SKIP_NOSOURCE=YES -wo NUM_THREADS=ALL_CPUS
> -wo OPTIMIZE_SIZE=TRUE -tr 1.0 1.0 -of vrt -s_srs ds_4326.prj -t_srs
> ds_4258.prj vrt_with_target_srs.vrt
> 
> That command generated this result:
> 
> GDAL:
> GDALOpen(/users/ds57201a/user3/rbda/de/uka30qs/dlu/grd/617bc9d0-52ce-11e3-
> 9de2-f1b38c98f22b/uka_parameter/vrt_with_source_srs.vrt, this=44f20)
> succeeds as VRT. OGRCT: PROJ >= 4.8.0 features enabled
> OGRCT: Wrap source at 6.16833.
> OGRCT: Source: +proj=longlat +ellps=WGS84 +no_defs
> OGRCT: Target: +proj=longlat +ellps=GRS80 +no_defs
> OGRCT: Wrap target at 6.16833.
> OGRCT: Source: +proj=longlat +ellps=GRS80 +no_defs
> OGRCT: Target: +proj=longlat +ellps=WGS84 +no_defs
> OGRCT: Wrap source at 6.16833.
> OGRCT: Source: +proj=longlat +ellps=WGS84 +no_defs
> OGRCT: Target: +proj=longlat +ellps=GRS80 +no_defs
> OGRCT: Wrap target at 6.16833.
> OGRCT: Source: +proj=longlat +ellps=GRS80 +no_defs
> OGRCT: Target: +proj=longlat +ellps=WGS84 +no_defs
> Creating output file that is 0P x 0L.
> ERROR 1: Attempt to create 0x0 dataset is illegal,sizes must be larger than
> zero.
> 
> Gdalinfo on the warped dataset(when I don't get the above error) shows
> this: Corner Coordinates:
> - Upper Left  (   30301.000,  187402.000) (Invalid angle,Invalid angle)
> Lower Left  (   30301.000,  178241.000) (Invalid angle,Invalid angle)
> Upper Right (   39346.000,  187402.000) (Invalid angle,Invalid angle)
> Lower Right (   39346.000,  178241.000) (Invalid angle,Invalid angle)
> Center      (   34823.500,  182821.500) (Invalid angle,Invalid angle)
> 
> Is there anyone who has experienced the same thing or  knows how to deal
> with data in decimal seconds within GDAL?
> 
> 
> Thanks and regards
> Moses

-- 
Geospatial professional services
http://even.rouault.free.fr/services.html


More information about the gdal-dev mailing list