Quantcast

[gdal-dev] GeoPDF translation

classic Classic list List threaded Threaded
13 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

[gdal-dev] GeoPDF translation

Smith, Michael-2
I am trying to convert the historic USGS topo maps (GeoPDFs) into GeoTIFFs  with GDAL 1.8 .  A very straightforward thing to do is
 
gdal_translate in.pdf out.tiff
 
This works fine but the TIFF looks as if the PDF was exported at 150dpi, which is a pretty crappy resolution.  Using something else (like PDFCreator and printing to a TIFF), the USGS GeoPDFs can be exported at higher resolutions, but then of course you lose the georeferencing.  So I have GDAL which keeps the georeferencing but only seems to like 150dpi, and then a zillion other solutions which allow 300 dpi or 600 dpi, but have no georeferencing.
 
I did try gdalwarp -tr 5 5 for example (which would roughly be equivalent to one of these maps at 300dpi), but it is clear that it is still rendering the PDF at 150dpi and then generating a 5-meter TIFF - looks no different really than the default output with gdal_translate.
 
Is there some argument that would specify the dpi at which the GeoPDF were rendered before conversion to GeoTIFF?  I don't see anything like that in the docs or the list archives.
 
BTW if you haven't checked out the historic topos they are great, going back to the 1880s for Maine.  What fun!
 
My ultimate goal is to have a seamless mosaic WMS of the old topos for Maine.
 
*****************
Michael Smith
State GIS Manager, Maine Office of GIS, Maine OIT
Board Member, Maine GeoLibrary
Board Member, Maine GIS Users Group
State Rep, National States Geographic Information Council

icons MEGISGeoLibraryMEGUGNSGIC

State House Station 174
264 Civic Center Drive
Augusta, ME 04333-0174
(207) 215-5530

69 o 47' 49.5"W   44 o 20' 54.5"N


 

_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev

icons_tiny.jpg (11K) Download Attachment
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: GeoPDF translation

Even Rouault
>
> Is there some argument that would specify the dpi at which the GeoPDF
> were rendered before conversion to GeoTIFF?  I don't see anything like
> that in the docs or the list archives.

Oh! Did you check http://gdal.org/frmt_pdf.html ? If so, check again ;-)

>
> BTW if you haven't checked out the historic topos they are great, going
> back to the 1880s for Maine.  What fun!
>
> My ultimate goal is to have a seamless mosaic WMS of the old topos for
> Maine.
>
> *****************
> Michael Smith
> State GIS Manager, Maine Office of GIS, Maine OIT
> Board Member, Maine GeoLibrary
> Board Member, Maine GIS Users Group
> State Rep, National States Geographic Information Council
>
>
>
> State House Station 174
> 264 Civic Center Drive
> Augusta, ME 04333-0174
> (207) 215-5530
>
> 69 o 47' 49.5"W   44 o 20' 54.5"N
_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: GeoPDF translation

Brent Fraser
In reply to this post by Smith, Michael-2
Michael,

  GDAL's PDF format page says you can use the GDAL_PDF_DPI config option...
Best Regards,
Brent Fraser

On 12/13/2011 11:58 AM, Smith, Michael wrote:
I am trying to convert the historic USGS topo maps (GeoPDFs) into GeoTIFFs  with GDAL 1.8 .  A very straightforward thing to do is
 
gdal_translate in.pdf out.tiff
 
This works fine but the TIFF looks as if the PDF was exported at 150dpi, which is a pretty crappy resolution.  Using something else (like PDFCreator and printing to a TIFF), the USGS GeoPDFs can be exported at higher resolutions, but then of course you lose the georeferencing.  So I have GDAL which keeps the georeferencing but only seems to like 150dpi, and then a zillion other solutions which allow 300 dpi or 600 dpi, but have no georeferencing.
 
I did try gdalwarp -tr 5 5 for example (which would roughly be equivalent to one of these maps at 300dpi), but it is clear that it is still rendering the PDF at 150dpi and then generating a 5-meter TIFF - looks no different really than the default output with gdal_translate.
 
Is there some argument that would specify the dpi at which the GeoPDF were rendered before conversion to GeoTIFF?  I don't see anything like that in the docs or the list archives.
 
BTW if you haven't checked out the historic topos they are great, going back to the 1880s for Maine.  What fun!
 
My ultimate goal is to have a seamless mosaic WMS of the old topos for Maine.
 
*****************
Michael Smith
State GIS Manager, Maine Office of GIS, Maine OIT
Board Member, Maine GeoLibrary
Board Member, Maine GIS Users Group
State Rep, National States Geographic Information Council

icons MEGIS GeoLibrary MEGUG NSGIC

State House Station 174
264 Civic Center Drive
Augusta, ME 04333-0174
(207) 215-5530

69 o 47' 49.5"W   44 o 20' 54.5"N


 


_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev

_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: GeoPDF translation

Derrick J Brashear
In reply to this post by Smith, Michael-2
>From my mail to geowankers:


gdalinfo $i|grep NEATLINE | awk -F= '{print "foo,WKT\nbla,\"" $2 "\""}' > $i.csv;
gdalwarp -crop_to_cutline -cutline $i.csv  -co "GDAL_PDF_DPI=250" -of GTiff  $i $i.tiff


On Tue, Dec 13, 2011 at 1:58 PM, Smith, Michael <[hidden email]> wrote:
I am trying to convert the historic USGS topo maps (GeoPDFs) into GeoTIFFs  with GDAL 1.8 .  A very straightforward thing to do is
 
gdal_translate in.pdf out.tiff
 
This works fine but the TIFF looks as if the PDF was exported at 150dpi, which is a pretty crappy resolution.  Using something else (like PDFCreator and printing to a TIFF), the USGS GeoPDFs can be exported at higher resolutions, but then of course you lose the georeferencing.  So I have GDAL which keeps the georeferencing but only seems to like 150dpi, and then a zillion other solutions which allow 300 dpi or 600 dpi, but have no georeferencing.
 
I did try gdalwarp -tr 5 5 for example (which would roughly be equivalent to one of these maps at 300dpi), but it is clear that it is still rendering the PDF at 150dpi and then generating a 5-meter TIFF - looks no different really than the default output with gdal_translate.
 
Is there some argument that would specify the dpi at which the GeoPDF were rendered before conversion to GeoTIFF?  I don't see anything like that in the docs or the list archives.
 
BTW if you haven't checked out the historic topos they are great, going back to the 1880s for Maine.  What fun!
 
My ultimate goal is to have a seamless mosaic WMS of the old topos for Maine.
 
*****************
Michael Smith
State GIS Manager, Maine Office of GIS, Maine OIT
Board Member, Maine GeoLibrary
Board Member, Maine GIS Users Group
State Rep, National States Geographic Information Council

icons MEGISGeoLibraryMEGUGNSGIC

State House Station 174
264 Civic Center Drive
Augusta, ME 04333-0174
<a href="tel:%28207%29%20215-5530" value="+12072155530" target="_blank">(207) 215-5530

69 o 47' 49.5"W   44 o 20' 54.5"N


 

_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev



--
Derrick

_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: GeoPDF translation

Derrick J Brashear


On Tue, Dec 13, 2011 at 2:13 PM, Derrick Brashear <[hidden email]> wrote:
From my mail to geowankers:


gdalinfo $i|grep NEATLINE | awk -F= '{print "foo,WKT\nbla,\"" $2 "\""}' > $i.csv;
gdalwarp -crop_to_cutline -cutline $i.csv  -co "GDAL_PDF_DPI=250" -of GTiff  $i $i.tiff


On Tue, Dec 13, 2011 at 1:58 PM, Smith, Michael <[hidden email]> wrote:
I am trying to convert the historic USGS topo maps (GeoPDFs) into GeoTIFFs  with GDAL 1.8 .  A very straightforward thing to do is
 

Oh, and you want 1.8.1 for -crop_to_cutline

none of this is needed if you're not tiling; this is just the extras to crop to the neatline. otherwise, the GDAL_PDF_DPI=250 (which is the native scanning resolution per the usgs)
is the relevant bit.
 
gdal_translate in.pdf out.tiff
 
This works fine but the TIFF looks as if the PDF was exported at 150dpi, which is a pretty crappy resolution.  Using something else (like PDFCreator and printing to a TIFF), the USGS GeoPDFs can be exported at higher resolutions, but then of course you lose the georeferencing.  So I have GDAL which keeps the georeferencing but only seems to like 150dpi, and then a zillion other solutions which allow 300 dpi or 600 dpi, but have no georeferencing.
 
I did try gdalwarp -tr 5 5 for example (which would roughly be equivalent to one of these maps at 300dpi), but it is clear that it is still rendering the PDF at 150dpi and then generating a 5-meter TIFF - looks no different really than the default output with gdal_translate.
 
Is there some argument that would specify the dpi at which the GeoPDF were rendered before conversion to GeoTIFF?  I don't see anything like that in the docs or the list archives.
 
BTW if you haven't checked out the historic topos they are great, going back to the 1880s for Maine.  What fun!
 
My ultimate goal is to have a seamless mosaic WMS of the old topos for Maine.
 
*****************
Michael Smith
State GIS Manager, Maine Office of GIS, Maine OIT
Board Member, Maine GeoLibrary
Board Member, Maine GIS Users Group
State Rep, National States Geographic Information Council

icons MEGISGeoLibraryMEGUGNSGIC

State House Station 174
264 Civic Center Drive
Augusta, ME 04333-0174
<a href="tel:%28207%29%20215-5530" value="+12072155530" target="_blank">(207) 215-5530

69 o 47' 49.5"W   44 o 20' 54.5"N


 

_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev



--
Derrick



--
Derrick

_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

RE: GeoPDF translation

Smith, Michael-2
Thanks to those of you who pointed out the obvious -co "GDAL_PDF_DPI" which I somehow overlooked in the online docs.  RTFM right...
 
That said, I ran a test script using that CO and GDAL yields a warning message and then makes the export at 150dpi anyway.
 
gdal_translate -co "GDAL_PDF_DPI=300" -co "TFW=YES" addison.pdf addison.tif
 
This is returned:
Warning 6: Driver GTiff does not support GDAL_PDF_DPI creation option
 
and then GDAL goes ahead and translates the PDF to GTIFF using the default 150dpi.
 
*****************
Michael Smith
State GIS Manager, Maine Office of GIS, Maine OIT
Board Member, Maine GeoLibrary
Board Member, Maine GIS Users Group
State Rep, National States Geographic Information Council

icons MEGISGeoLibraryMEGUGNSGIC

State House Station 174
264 Civic Center Drive
Augusta, ME 04333-0174
(207) 215-5530

69 o 47' 49.5"W   44 o 20' 54.5"N


 

_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: GeoPDF translation

Even Rouault
Le mardi 13 décembre 2011 21:18:51, Smith, Michael a écrit :

> Thanks to those of you who pointed out the obvious -co "GDAL_PDF_DPI"
> which I somehow overlooked in the online docs.  RTFM right...
>
> That said, I ran a test script using that CO and GDAL yields a warning
> message and then makes the export at 150dpi anyway.
>
> gdal_translate -co "GDAL_PDF_DPI=300" -co "TFW=YES" addison.pdf
> addison.tif
>
> This is returned:
> Warning 6: Driver GTiff does not support GDAL_PDF_DPI creation option
>
> and then GDAL goes ahead and translates the PDF to GTIFF using the
> default 150dpi.

Yes,  GDAL_PDF_DPI is not a creation option of the GTiff driver, but a
configuration option of the PDF driver.

Configuration options are specified with "--config optionname optionvalue" (note
: space between name and value, not '=' character)

So :

gdal_translate --config GDAL_PDF_DPI 300 -co "TFW=YES" addison.pdf addison.tif
_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: GeoPDF translation

jsschmidt
This post has NOT been accepted by the mailing list yet.
Here's a Python3.3 script to convert multiple historic USGS GeoPDF's to GeoTiffs in a UTM map projection. (Script can be easily modified to output to other map projections). Intermediate Geotiffs are output in the native map projection - polyconic - then converted to GeoTiffs in UTM.  Script requires access to the GDAL utilities, which I obtained by installing MS4W and running the script from the MS4W command prompt.
Modify the script to point to the directory where your historic GeoPDFs are stored.

#Python 3.3 script to convert historic USGS PDF maps to GeoTiff format in UTM coordinates.
#Relies on GDAL (Geospatial Data Abstraction Library)
import glob, os
import subprocess
import pyproj

def GetLatLon(line):
    coords = line.split(') (')[1]
    coords = coords[:-1]
    LonStr, LatStr = coords.split(',')
    # Longitude
    LonStr = LonStr.split('d')    # Get the degrees, and the rest
    LonD = int(LonStr[0])
    LonStr = LonStr[1].split('\'')# Get the arc-m, and the rest
    LonM = int(LonStr[0])
    LonStr = LonStr[1].split('"') # Get the arc-s, and the rest
    LonS = float(LonStr[0])
    Lon = LonD + LonM/60. + LonS/3600.
    if LonStr[1] in ['W', 'w']:
        Lon = -1*Lon
    # Same for Latitude
    LatStr = LatStr.split('d')
    LatD = int(LatStr[0])
    LatStr = LatStr[1].split('\'')
    LatM = int(LatStr[0])
    LatStr = LatStr[1].split('"')
    LatS = float(LatStr[0])
    Lat = LatD + LatM/60. + LatS/3600.
    if LatStr[1] in ['S', 's']:
        Lat = -1*Lat
    return Lat, Lon

EPSGin = "EPSG:4267"                   # Lat-Lon, NAD27 - Assumed coordinate system of lat-lon corner coordinates as defined by EPSG code
EPSGout = "EPSG:26910"                 # UTM10, NAD83 - Output coordinate system for Geotiff files as defined by EPSG code
#EPSGout = "EPSG:26911"                 # UTM11, NAD83 - Output coordinate system for Geotiff files as defined by EPSG code

dirname = "C:\\GeoPDF\\"  # Location of input pdf files Also used for output directory
print(dirname)

p1=pyproj.Proj(init= EPSGin)           # Define input projection
p2=pyproj.Proj(init= EPSGout)          # Define output projection

list_of_files1 = glob.glob(dirname +'*.pdf')           # create the list of PDF files


for file_name in list_of_files1:
  os.rename(file_name, file_name.replace(" ", "_"))    #replace any blanks in PDF filenames with underscore character

list_of_files2 = glob.glob(dirname +'*.pdf')           # create the list of renamed PDF files

for file_name in list_of_files2:

  FIPDF = file_name #Input GeoPDF filename
  FOTIFF1 = FIPDF.replace('pdf', 'tif') # Output Filename for intermediate GeoTiff file in original map coordinates
  FOTIFF2 = FOTIFF1.replace('CA_', 'CA_UTM_') # Output Filename for final GeoTiff file in UTM Coordinates
  print(FIPDF)
  print(FOTIFF1)
  print(FOTIFF2)

# Export geoPDF to Geotiff at 300 dpi in original map projection coordinates (polyconic for USGS historic maps)
# Use JPEG Compression in the GeoTiff

  os.system('gdal_translate.exe  -of GTiff --config GDAL_PDF_DPI 300 -co COMPRESS=JPEG ' + FIPDF + ' '+ FOTIFF1)

  CornerLats = [0,0,0,0]
  CornerLons = [0,0,0,0]
  Cornerx = [0,0,0,0]
  Cornery = [0,0,0,0]
  gcp = ['','','','']
  xstr = ['','','','']
  ystr = ['','','','']

# Run GDALinfo to read the header information from each intermediate GeoTiff file

  GdalInfo = subprocess.check_output('gdalinfo {}'.format(FOTIFF1),shell=True,universal_newlines=True)# Runs Gdalinfo on exported GeoTiff
  GdalInfo = GdalInfo.splitlines() # Creates a line by line list.
  GotUL, GotUR, GotLL, GotLR, GotSize, GotPixel, Proj, NeedReproj = False, False, False, False, False, False, False, True
 
  for line in GdalInfo:
        print(line)  
        if line.startswith("Size"):
# Get the row-column size of the GeoTiff
            size = line.split()
            R = size[3]
            size2 = size[2].split(',')
            C = size2[0]
            CornerRow = ['0', R, '0', R]
            CornerCol = ['0', '0', C, C]
            GotSize = True
       
        if line.startswith("Pixel Size"):
# Calculate the size of each pixel in map units (assumed to be meters)
            pix1 = line.split('(')
            pix2 = pix1[1].split(',')
            pixsize = float(pix2[0])
            halfpix = pixsize/2 # Adjustment to move coordinates to cell center from edge
            GotPixel = True
             
        if line.startswith("Upper Left"):
            CornerLats[0], CornerLons[0] = GetLatLon(line)
            Cornerx[0], Cornery[0] = pyproj.transform(p1,p2,CornerLons[0],CornerLats[0]) #Reproject edge coords in UL to output Map proj
            Cornerx[0] = Cornerx[0] - halfpix # Adjustment to pixel edge
            Cornery[0] = Cornery[0] + halfpix
            GotUL = True
       
        if line.startswith("Lower Left"):
            CornerLats[1], CornerLons[1] = GetLatLon(line)
            Cornerx[1], Cornery[1] = pyproj.transform(p1,p2,CornerLons[1],CornerLats[1]) #Reproject edge coords in LL to output Map proj
            Cornerx[1] = Cornerx[1] - halfpix
            Cornery[1] = Cornery[1] - halfpix
            GotLL = True

        if line.startswith("Upper Right"):
            CornerLats[2], CornerLons[2] = GetLatLon(line)
            Cornerx[2], Cornery[2] = pyproj.transform(p1,p2,CornerLons[2],CornerLats[2]) #Reproject edge coords in UR to output Map proj
            Cornerx[2] = Cornerx[2] + halfpix
            Cornery[2] = Cornery[2] + halfpix
            GotUR = True

        if line.startswith("Lower Right"):
            CornerLats[3], CornerLons[3] = GetLatLon(line)
            Cornerx[3], Cornery[3] = pyproj.transform(p1,p2,CornerLons[3],CornerLats[3]) #Reproject edge coords in LR to output Map proj
            Cornerx[3] = Cornerx[3] + halfpix
            Cornery[3] = Cornery[3] - halfpix
            GotLR = True

        if GotUL and GotUR and GotLL and GotLR and GotSize and GotPixel and NeedReproj:
            xstr[0] = str(Cornerx[0])
            xstr[1] = str(Cornerx[1])
            xstr[2] = str(Cornerx[2])
            xstr[3] = str(Cornerx[3])
            ystr[0] = str(Cornery[0])
            ystr[1] = str(Cornery[1])
            ystr[2] = str(Cornery[2])
            ystr[3] = str(Cornery[3])

            gcp[0] = ' -gcp '+ CornerCol[0]+ ' ' + CornerRow[0]+ ' ' + xstr[0]+ ' ' + ystr[0]
            gcp[1] = ' -gcp '+ CornerCol[1]+ ' ' + CornerRow[1]+ ' ' + xstr[1]+ ' ' + ystr[1]
            gcp[2] = ' -gcp '+ CornerCol[2]+ ' ' + CornerRow[2]+ ' ' + xstr[2]+ ' ' + ystr[2]
            gcp[3] = ' -gcp '+ CornerCol[3]+ ' ' + CornerRow[3]+ ' ' + xstr[3]+ ' ' + ystr[3]

            # Reproject to new geotiff file using control points
 
            gdalstr = 'gdal_translate.exe -of GTIFF -co COMPRESS=JPEG ' + FOTIFF1 + gcp[0] + gcp[1] + gcp[2] + gcp[3]+ ' -a_srs ' + EPSGout + ' '+ FOTIFF2
            print(gdalstr)
            os.system(gdalstr)
            NeedReproj = False
            Proj = True

        if GotUL and GotUR and GotLL and GotLR and GotSize and GotPixel and Proj:
            break

         
           
 



Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: [gdal-dev] GeoPDF translation

jsschmidt
In reply to this post by Smith, Michael-2

Here's a Python3.3 script to convert multiple historic USGS GeoPDF's to GeoTiffs in a UTM10, NAD83 map projection at 300 dpi. Intermediate Geotiffs are output in the native USGS map projection - polyconic - then converted to GeoTiffs in UTM.  Script requires access to the GDAL utilities, which I obtained by installing MS4W and running the script from the MS4W command prompt.
Modify the script to point to the directory where your historic GeoPDFs are stored, to change the output map projection, and to modify the name of the output GeoTiff files. (Script currently assumes PDF files begin with 'CA' and final output files will be named CA_UTM....tif.)









#Python 3.3 script to convert historic USGS PDF maps to GeoTiff format in UTM coordinates.
#Relies on GDAL (Geospatial Data Abstraction Library)
import glob, os
import subprocess
import pyproj

def GetLatLon(line):
    coords = line.split(') (')[1]
    coords = coords[:-1]
    LonStr, LatStr = coords.split(',')
    # Longitude
    LonStr = LonStr.split('d')    # Get the degrees, and the rest
    LonD = int(LonStr[0])
    LonStr = LonStr[1].split('\'')# Get the arc-m, and the rest
    LonM = int(LonStr[0])
    LonStr = LonStr[1].split('"') # Get the arc-s, and the rest
    LonS = float(LonStr[0])
    Lon = LonD + LonM/60. + LonS/3600.
    if LonStr[1] in ['W', 'w']:
        Lon = -1*Lon
    # Same for Latitude
    LatStr = LatStr.split('d')
    LatD = int(LatStr[0])
    LatStr = LatStr[1].split('\'')
    LatM = int(LatStr[0])
    LatStr = LatStr[1].split('"')
    LatS = float(LatStr[0])
    Lat = LatD + LatM/60. + LatS/3600.
    if LatStr[1] in ['S', 's']:
        Lat = -1*Lat
    return Lat, Lon

EPSGin = "EPSG:4267"                   # Lat-Lon, NAD27 - Assumed coordinate system of lat-lon corner coordinates as defined by EPSG code
EPSGout = "EPSG:26910"                 # UTM10, NAD83 - Output coordinate system for Geotiff files as defined by EPSG code
#EPSGout = "EPSG:26911"                 # UTM11, NAD83 - Output coordinate system for Geotiff files as defined by EPSG code

dirname = "C:\\GeoPDF\\"  # Location of input pdf files Also used for output directory
print(dirname)

p1=pyproj.Proj(init= EPSGin)           # Define input projection
p2=pyproj.Proj(init= EPSGout)          # Define output projection

list_of_files1 = glob.glob(dirname +'*.pdf')           # create the list of PDF files


for file_name in list_of_files1:
  os.rename(file_name, file_name.replace(" ", "_"))    #replace any blanks in PDF filenames with underscore character

list_of_files2 = glob.glob(dirname +'*.pdf')           # create the list of renamed PDF files

for file_name in list_of_files2:

  FIPDF = file_name #Input GeoPDF filename
  FOTIFF1 = FIPDF.replace('pdf', 'tif') # Output Filename for intermediate GeoTiff file in original map coordinates
  FOTIFF2 = FOTIFF1.replace('CA_', 'CA_UTM_') # Output Filename for final GeoTiff file in UTM Coordinates
  print(FIPDF)
  print(FOTIFF1)
  print(FOTIFF2)

# Export geoPDF to Geotiff at 300 dpi in original map projection coordinates (polyconic for USGS historic maps)
# Use JPEG Compression in the GeoTiff

  os.system('gdal_translate.exe  -of GTiff --config GDAL_PDF_DPI 300 -co COMPRESS=JPEG ' + FIPDF + ' '+ FOTIFF1)

  CornerLats = [0,0,0,0]
  CornerLons = [0,0,0,0]
  Cornerx = [0,0,0,0]
  Cornery = [0,0,0,0]
  gcp = ['','','','']
  xstr = ['','','','']
  ystr = ['','','','']

# Run GDALinfo to read the header information from each intermediate GeoTiff file

  GdalInfo = subprocess.check_output('gdalinfo {}'.format(FOTIFF1),shell=True,universal_newlines=True)# Runs Gdalinfo on exported GeoTiff
  GdalInfo = GdalInfo.splitlines() # Creates a line by line list.
  GotUL, GotUR, GotLL, GotLR, GotSize, GotPixel, Proj, NeedReproj = False, False, False, False, False, False, False, True
 
  for line in GdalInfo:
        print(line)  
        if line.startswith("Size"):
# Get the row-column size of the GeoTiff
            size = line.split()
            R = size[3]
            size2 = size[2].split(',')
            C = size2[0]
            CornerRow = ['0', R, '0', R]
            CornerCol = ['0', '0', C, C]
            GotSize = True
       
        if line.startswith("Pixel Size"):
# Calculate the size of each pixel in map units (assumed to be meters)
            pix1 = line.split('(')
            pix2 = pix1[1].split(',')
            pixsize = float(pix2[0])
            halfpix = pixsize/2 # Adjustment to move coordinates to cell center from edge
            GotPixel = True
             
        if line.startswith("Upper Left"):
            CornerLats[0], CornerLons[0] = GetLatLon(line)
            Cornerx[0], Cornery[0] = pyproj.transform(p1,p2,CornerLons[0],CornerLats[0]) #Reproject edge coords in UL to output Map proj
            Cornerx[0] = Cornerx[0] - halfpix # Adjustment to pixel edge
            Cornery[0] = Cornery[0] + halfpix
            GotUL = True
       
        if line.startswith("Lower Left"):
            CornerLats[1], CornerLons[1] = GetLatLon(line)
            Cornerx[1], Cornery[1] = pyproj.transform(p1,p2,CornerLons[1],CornerLats[1]) #Reproject edge coords in LL to output Map proj
            Cornerx[1] = Cornerx[1] - halfpix
            Cornery[1] = Cornery[1] - halfpix
            GotLL = True

        if line.startswith("Upper Right"):
            CornerLats[2], CornerLons[2] = GetLatLon(line)
            Cornerx[2], Cornery[2] = pyproj.transform(p1,p2,CornerLons[2],CornerLats[2]) #Reproject edge coords in UR to output Map proj
            Cornerx[2] = Cornerx[2] + halfpix
            Cornery[2] = Cornery[2] + halfpix
            GotUR = True

        if line.startswith("Lower Right"):
            CornerLats[3], CornerLons[3] = GetLatLon(line)
            Cornerx[3], Cornery[3] = pyproj.transform(p1,p2,CornerLons[3],CornerLats[3]) #Reproject edge coords in LR to output Map proj
            Cornerx[3] = Cornerx[3] + halfpix
            Cornery[3] = Cornery[3] - halfpix
            GotLR = True

        if GotUL and GotUR and GotLL and GotLR and GotSize and GotPixel and NeedReproj:
            xstr[0] = str(Cornerx[0])
            xstr[1] = str(Cornerx[1])
            xstr[2] = str(Cornerx[2])
            xstr[3] = str(Cornerx[3])
            ystr[0] = str(Cornery[0])
            ystr[1] = str(Cornery[1])
            ystr[2] = str(Cornery[2])
            ystr[3] = str(Cornery[3])

            gcp[0] = ' -gcp '+ CornerCol[0]+ ' ' + CornerRow[0]+ ' ' + xstr[0]+ ' ' + ystr[0]
            gcp[1] = ' -gcp '+ CornerCol[1]+ ' ' + CornerRow[1]+ ' ' + xstr[1]+ ' ' + ystr[1]
            gcp[2] = ' -gcp '+ CornerCol[2]+ ' ' + CornerRow[2]+ ' ' + xstr[2]+ ' ' + ystr[2]
            gcp[3] = ' -gcp '+ CornerCol[3]+ ' ' + CornerRow[3]+ ' ' + xstr[3]+ ' ' + ystr[3]

            # Reproject to new geotiff file using control points
 
            gdalstr = 'gdal_translate.exe -of GTIFF -co COMPRESS=JPEG ' + FOTIFF1 + gcp[0] + gcp[1] + gcp[2] + gcp[3]+ ' -a_srs ' + EPSGout + ' '+ FOTIFF2
            print(gdalstr)
            os.system(gdalstr)
            NeedReproj = False
            Proj = True

        if GotUL and GotUR and GotLL and GotLR and GotSize and GotPixel and Proj:
            break

         
           
 



Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

[gdal-dev] Passing RAT between drivers, including VRT

Ivan Lucena-4
Hi folks,

I might be wrong, but it looks like the VRT driver is not supporting RAT at this moment (TRUNK).

Is that correct?

A simple way to test it is by "gdal_translating" that file "i8u_c_i.img" from autotest to VRT:

$ gdal_translate -of vrt i8u_c_i.img vrt_i8u_c_i.vrt
$ gdalinfo vrt_i8u_c_i.vrt -noct

The RAT is not there.

That wouldn't be a big problem by itself, the real problem is that gdal_translate uses an internal VRT to support several options. And I am trying to fix that ticket: http://trac.osgeo.org/gdal/ticket/5200 and I noticed that when I use *-a_nodata* the GetDefaultRAT() it returns null.

Drivers that implement *file based* GDALPamDataset doesn't seems to have that problem because they don't need to call GetDefaultRAT(). Examples:

$ gdal_translate -of gtiif i8u_c_i.img gtiff_i8u_c_i.tif -a_nodata 0
$ gdal_translate -of gtiif i8u_c_i.img gtiff_i8u_c_i.tif

The RAT is transferred to the geotiff file in both instances.

Is that something we would like to get fix in VRT?

Best regards,

Ivan




_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Passing RAT between drivers, including VRT

Even Rouault
Le jeudi 15 août 2013 18:53:46, Ivan Lucena a écrit :
> Hi folks,
>
> I might be wrong, but it looks like the VRT driver is not supporting RAT at
> this moment (TRUNK).
>
> Is that correct?

Hi Ivan,

yes, correct. See http://trac.osgeo.org/gdal/ticket/4903 where there's a patch
pending for this support.

>
> A simple way to test it is by "gdal_translating" that file "i8u_c_i.img"
> from autotest to VRT:
>
> $ gdal_translate -of vrt i8u_c_i.img vrt_i8u_c_i.vrt
> $ gdalinfo vrt_i8u_c_i.vrt -noct
>
> The RAT is not there.
>
> That wouldn't be a big problem by itself, the real problem is that
> gdal_translate uses an internal VRT to support several options. And I am
> trying to fix that ticket: http://trac.osgeo.org/gdal/ticket/5200 and I
> noticed that when I use *-a_nodata* the GetDefaultRAT() it returns null.
>
> Drivers that implement *file based* GDALPamDataset doesn't seems to have
> that problem because they don't need to call GetDefaultRAT(). Examples:
>
> $ gdal_translate -of gtiif i8u_c_i.img gtiff_i8u_c_i.tif -a_nodata 0
> $ gdal_translate -of gtiif i8u_c_i.img gtiff_i8u_c_i.tif
>
> The RAT is transferred to the geotiff file in both instances.

That's surprising !, and even not what I observe in my test. The RAT is
transfered in the simple gdal_translate case (since the IMG dataset is
directly passed as the source dataset so GetDefaultRAT() != NULL), but not in
the -a_nodata 0 case since the source dataset is a VRT one, without RAT
support.

>
> Is that something we would like to get fix in VRT?

Likely, the patches in http://trac.osgeo.org/gdal/ticket/4903 would need to be
rebased on trunk since RFC 40 has been implemented.

>
> Best regards,
>
> Ivan

--
Geospatial professional services
http://even.rouault.free.fr/services.html
_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Passing RAT between drivers, including VRT

Ivan Lucena-4
Hi Even,

I double checked and I noticed that this is not true:

> > The RAT is transferred to the geotiff file in both instances.

I pass -a_nodata and got no RAT as you can see on the following log.

Without -a_nodata the file gtiff_i8u_c_i.tif.aux.xml is create with the full RAT.

$ gdal_translate -of gtiff i8u_c_i.img gtiff_i8u_c_i.tif -a_nodata 0
Input file size is 233, 250
0...10...20...30...40...50...60...70...80...90...100 - done.

$ gdalinfo gtiff_i8u_c_i.tif -noct
Driver: GTiff/GeoTIFF
Files: gtiff_i8u_c_i.tif
Size is 233, 250
Coordinate System is `'
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  250.0)
Upper Right (  233.0,    0.0)
Lower Right (  233.0,  250.0)
Center      (  116.5,  125.0)
Band 1 Block=233x35 Type=Byte, ColorInterp=Palette
  Description = Band_1
  Min=0.000 Max=255.000
  Minimum=0.000, Maximum=255.000, Mean=118.537, StdDev=102.769
  NoData Value=0
  Metadata:
    LAYER_TYPE=thematic
    STATISTICS_HISTOBINVALUES=12603|1|0|0|45|1|0|0|0|0|656|177|0|0|5026|1062|0|0
|2|0|0|0|0|0|0|0|0|0|0|0|0|0|75|1|0|0|207|158|0|0|8|34|0|0|0|0|538|57|0|10|214|2
0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|1|31|0|0|9|625|67|0|0|118|738|117|3004|1499|49
1|187|1272|513|1|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|16|3|0|0|283|123|5|1931|835
|357|332|944|451|80|40|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|12|5|0|0|535|10
29|118|0|33|246|342|0|0|10|8|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|169|439|0
|0|6|990|329|0|0|120|295|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|164
|42|0|0|570|966|0|0|18|152|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|4
5|106|0|0|16|16517|
    STATISTICS_HISTOMAX=255
    STATISTICS_HISTOMIN=0
    STATISTICS_HISTONUMBINS=256
    STATISTICS_MAXIMUM=255
    STATISTICS_MEAN=118.53732188841
    STATISTICS_MEDIAN=85
    STATISTICS_MINIMUM=0
    STATISTICS_MODE=255
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_STDDEV=102.76909129768
  Color Table (RGB with 256 entries)

But for output HFA that is still true. I can pass -a_nodata or not and the RAT is always copied to the output.

Try that:

$ gdal_translate -of hfa i8u_c_i.img hfa_i8u_c_i.img -a_nodata 0
$ gdal_translate -of hfa i8u_c_i.img hfa_i8u_c_i.img

Is that the same that you are getting over there?

Regards,

Ivan


> From: [hidden email]

> To: [hidden email]
> Subject: Re: [gdal-dev] Passing RAT between drivers, including VRT
> Date: Thu, 15 Aug 2013 19:28:51 +0200
> CC: [hidden email]; [hidden email]
>
> Le jeudi 15 août 2013 18:53:46, Ivan Lucena a écrit :
> > Hi folks,
> >
> > I might be wrong, but it looks like the VRT driver is not supporting RAT at
> > this moment (TRUNK).
> >
> > Is that correct?
>
> Hi Ivan,
>
> yes, correct. See http://trac.osgeo.org/gdal/ticket/4903 where there's a patch
> pending for this support.
>
> >
> > A simple way to test it is by "gdal_translating" that file "i8u_c_i.img"
> > from autotest to VRT:
> >
> > $ gdal_translate -of vrt i8u_c_i.img vrt_i8u_c_i.vrt
> > $ gdalinfo vrt_i8u_c_i.vrt -noct
> >
> > The RAT is not there.
> >
> > That wouldn't be a big problem by itself, the real problem is that
> > gdal_translate uses an internal VRT to support several options. And I am
> > trying to fix that ticket: http://trac.osgeo.org/gdal/ticket/5200 and I
> > noticed that when I use *-a_nodata* the GetDefaultRAT() it returns null.
> >
> > Drivers that implement *file based* GDALPamDataset doesn't seems to have
> > that problem because they don't need to call GetDefaultRAT(). Examples:
> >
> > $ gdal_translate -of gtiif i8u_c_i.img gtiff_i8u_c_i.tif -a_nodata 0
> > $ gdal_translate -of gtiif i8u_c_i.img gtiff_i8u_c_i.tif
> >
> > The RAT is transferred to the geotiff file in both instances.
>
> That's surprising !, and even not what I observe in my test. The RAT is
> transfered in the simple gdal_translate case (since the IMG dataset is
> directly passed as the source dataset so GetDefaultRAT() != NULL), but not in
> the -a_nodata 0 case since the source dataset is a VRT one, without RAT
> support.
>
> >
> > Is that something we would like to get fix in VRT?
>
> Likely, the patches in http://trac.osgeo.org/gdal/ticket/4903 would need to be
> rebased on trunk since RFC 40 has been implemented.
>
> >
> > Best regards,
> >
> > Ivan
>
> --
> Geospatial professional services
> http://even.rouault.free.fr/services.html

_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Passing RAT between drivers, including VRT

Even Rouault
Le jeudi 15 août 2013 20:04:51, Ivan Lucena a écrit :

> Hi Even,
>
> I double checked and I noticed that this is not true:
> > > The RAT is transferred to the geotiff file in both instances.
>
> I pass -a_nodata and got no RAT as you can see on the
> following log.
>
> Without -a_nodata the file gtiff_i8u_c_i.tif.aux.xml is create with the
> full RAT.

Yes, that's exactly what I said in my previous email.

>
> But for output HFA that is still true. I can pass
>  -a_nodata or not and the RAT is always copied to the output.
>
> Try that:
>
> $ gdal_translate -of hfa i8u_c_i.img hfa_i8u_c_i.img -a_nodata 0
> $ gdal_translate -of hfa i8u_c_i.img hfa_i8u_c_i.img
>
> Is that the same that you are getting over there?

Well, this is different. The HFA driver is a bit particular in its behaviour
and reports a RAT when there's a color table (in fact a color table is a
particular type of RAT). So technically the RAT of the source IMG is not
transferred into the target IMG in the -a_nodata 0 case, but the color table
is.


--
Geospatial professional services
http://even.rouault.free.fr/services.html
_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Loading...