Dear all,
Overall, I'm trying to convert ESRI SHAPE data, given in "Web Mercator"
EPSG:3857 , to a TIF in ETRS89 / UTM Zone 32U , EPSG:25832 .
The GDAL Version is the one included in current stable Debian,
https://packages.debian.org/bookworm/gdal-bin , gdal-bin (3.6.2+dfsg-1 and others)
The first step, rasterization, works flawless:
gdal_rasterize -a_srs EPSG:3857 -tr 10 10 -burn 1 q.shp q.tif
gdalinfo lists the file as: (full text attached)
....
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
BASEGEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
..
Corner Coordinates:
Upper Left ( 836370.121, 6080896.490) ( 7d30'47.67"E, 47d50'37.43"N)
Lower Left ( 836370.121, 6045996.490) ( 7d30'47.67"E, 47d37'58.40"N)
Upper Right ( 887980.121, 6080896.490) ( 7d58'36.70"E, 47d50'37.43"N)
Lower Right ( 887980.121, 6045996.490) ( 7d58'36.70"E, 47d37'58.40"N)
Center ( 862175.121, 6063446.490) ( 7d44'42.18"E, 47d44'18.30"N)
Band 1 Block=5161x1 Type=Float64, ColorInterp=Gray
Running the lower left corner through gdaltransform gives plausible ETRS89
coordinates :
echo 836375.120800 6045997.312000 | gdaltransform -s_srs EPSG:3857 -t_srs
EPSG:25832
388314.422838145 5276569.91979615 0
However, running the TIF through gdalwarp, and checking the result with
gdalinfo, does not show a plausible bounding box:
gdalwarp -s_srs EPSG:3857 -t_srs EPSG:25832 q.tif w.tif
gdalinfo w.tif
...
Corner Coordinates:
Upper Left ( 836370.121, 6080896.490) (Invalid angle,Invalid angle)
Lower Left ( 836370.121, 6045996.490) (Invalid angle,Invalid angle)
Upper Right ( 887980.121, 6080896.490) (Invalid angle,Invalid angle)
Lower Right ( 887980.121, 6045996.490) (Invalid angle,Invalid angle)
Center ( 862175.121, 6063446.490) (Invalid angle,Invalid angle)
Looks like it passes gdalwarp unchanged.
Apparently, I'm overlooking something obvious.
Any pointers much appreciated,
thanks, best regards, have a nice Sunday and thanks for GDAL !
Peter
Driver: GTiff/GeoTIFF
Files: q.tif
Size is 5161, 3490
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
BASEGEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["Popular Visualisation Pseudo-Mercator",
METHOD["Popular Visualisation Pseudo Mercator",
ID["EPSG",1024]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",0,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["easting (X)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["northing (Y)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Web mapping and visualisation."],
AREA["World between 85.06°S and 85.06°N."],
BBOX[-85.06,-180,85.06,180]],
ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (836370.120799999684095,6080896.489900000393391)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 836370.121, 6080896.490) ( 7d30'47.67"E, 47d50'37.43"N)
Lower Left ( 836370.121, 6045996.490) ( 7d30'47.67"E, 47d37'58.40"N)
Upper Right ( 887980.121, 6080896.490) ( 7d58'36.70"E, 47d50'37.43"N)
Lower Right ( 887980.121, 6045996.490) ( 7d58'36.70"E, 47d37'58.40"N)
Center ( 862175.121, 6063446.490) ( 7d44'42.18"E, 47d44'18.30"N)
Band 1 Block=5161x1 Type=Float64, ColorInterp=Gray
Driver: GTiff/GeoTIFF
Files: w.tif
Size is 5161, 3490
Coordinate System is:
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (836370.120799999684095,6080896.489900000393391)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 836370.121, 6080896.490) (Invalid angle,Invalid angle)
Lower Left ( 836370.121, 6045996.490) (Invalid angle,Invalid angle)
Upper Right ( 887980.121, 6080896.490) (Invalid angle,Invalid angle)
Lower Right ( 887980.121, 6045996.490) (Invalid angle,Invalid angle)
Center ( 862175.121, 6063446.490) (Invalid angle,Invalid angle)
Band 1 Block=5161x1 Type=Float64, ColorInterp=Gray
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
