4

I run the code below, that I would hope to return the parameters from projection file (.prj). However it doesn't work as I would expected. I realized that the it isn't working because the some features of .prj file. I tried with another .prj file and the code has worked well. What is the problem to my code?

import ogr

# use Shapefile driver
driver = ogr.GetDriverByName("ESRI Shapefile")
# reference Shapefile
shp = "/home/dogosousa/Downloads/lc8_ba.shp"
# open the file
ds = driver.Open(shp, 0)
# reference the only layer in a Shapefile
lyr = ds.GetLayer(0)

print(lyr.GetSpatialRef())

# projected coordinate system
proj_string = lyr.GetSpatialRef().GetAttrValue("PROJCS", 0)
# geographic coordinate system
geog_string = lyr.GetSpatialRef().GetAttrValue("GEOGCS", 0)
# EPSG Code if available
epsg = lyr.GetSpatialRef().GetAttrValue("AUTHORITY", 1)
# datum
datum = lyr.GetSpatialRef().GetAttrValue("DATUM", 0)

print(proj_string, geog_string, epsg,datum)

The .prj file that doesn't work with the code:

PROJCS["PROJCS["WGS 84 / UTM zone 24N",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["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-39],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32624"]]",UNIT["Meter",1]]

The .prj file that works with the code:

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]
6

The WKT string in that .prj is corrupt. I don't know what wrote it, but the PROJCS and UNITS elements are repeated.

i.e. (indicated by --> <--)

-->PROJCS["<--PROJCS["WGS 84 / UTM zone 24N",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["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-39],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32624"]]",-->UNIT["Meter",1]<--]

Removing the dupes works

PROJCS["WGS 84 / UTM zone 24N",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["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-39],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32624"]]"]

e.g.

from osgeo import osr

srs = osr.SpatialReference()

prj='PROJCS["WGS 84 / UTM zone 24N",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["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-39],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32624"]]"'

srs.ImportFromWkt(prj)

# projected coordinate system
proj_string = srs.GetAttrValue("PROJCS", 0)
# geographic coordinate system
geog_string = srs.GetAttrValue("GEOGCS", 0)
# EPSG Code if available
epsg = srs.GetAttrValue("AUTHORITY", 1)
# datum
datum = srs.GetAttrValue("DATUM", 0)

print(proj_string, geog_string, epsg,datum)

WGS 84 / UTM zone 24N WGS 84 32624 WGS_1984

You can catch these sort of issues in future by enabling gdal exceptions, i.e. gdal.UseExceptions(), ogr.UseExceptions(), osr.UseExceptions(), e.g:

osr.UseExceptions()

prj='PROJCS["PROJCS["WGS 84 / UTM zone 24N",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["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJE
CTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-39],PARAMETER["scale_factor",0.9996],PARAMETER["false_eastin
g",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32624"]]"
,UNIT["Meter",1]]'

s.ImportFromWkt(prj)

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
...
RuntimeError: OGR Error: Corrupt data
  • So. The lib that wrote the file is gdal-segment. The problem comes from it. I don't know why. I will report. – Diogo Caribé Apr 21 '18 at 23:07
  • I don't know how, it just writes the exact SRS it reads from the source dataset github.com/cbalint13/gdal-segment/blob/master/src/io/…. Might be your input raster or a gdal bug. – user2856 Apr 21 '18 at 23:32
  • My input data are a composition from landsat 8. I only compose the bands that I downloaded from USGS with gdal_merge.py github.com/dogosousa/preprocess2ta/blob/master/ComposeBands.py. I cliped them too in order to test more fast. I think it wouldn't a problem. – Diogo Caribé Apr 22 '18 at 0:05
  • You can run gdalsrsinfo on any file in your workflow to see where GDAL refuses to read the CRS correctly. Might be a bug worth reporting if it is a GDAL component. – AndreJ Apr 22 '18 at 8:34
  • @AndreJ i will see it. – Diogo Caribé Apr 26 '18 at 0:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.