OGR_t1

OGR

Commands

  • ogrinfo: Lists information about an OGR-supported data source.
  • ogr2ogr: Converts simple features data between file formats.
  • ogrtindex: Creates a tileindex.
  • ogrlineref: Create linear reference and provide some calculations using it.
  • ogrmerge.py: Merge several vector datasets into a single one.
  • ogr_layer_algebra.py: Performs various Vector layer algebraic operations.
In [5]:
cd /Users/abharathi/Documents/gis_data/gdal-tools
/Users/abharathi/Documents/gis_data/gdal-tools
In [6]:
ls
1870_southern-india.jpg*           india-with-gcp.vrt
1870_southern-india.jpg.aux.xml    landsat8/
5d16a93f1cf0f6000579ad2c.tif       london_1m_dsm/
batch.py*                          naip/
batch_parallel.py*                 precipitation.gpkg
earth_at_night.jpg*                prism/
earth_at_night.tif                 reservoirs.xml
earthquakes/                       reservoirs_india.tif
gdal_stuff.qgz                     reservoirs_india_proj.tif
geonames/                          spatial_query.gpkg
india-reprojected-polynomial.tif   srtm/
india-reprojected-tps.tif          worldcities.csv*
india-reprojected-tps.tif.aux.xml

Info using OGRINFO

In [7]:
!ogrinfo worldcities.csv
INFO: Open of `worldcities.csv'
      using driver `CSV' successful.
1: worldcities (None)
In [9]:
!ogrinfo -al -so worldcities.csv
INFO: Open of `worldcities.csv'
      using driver `CSV' successful.

Layer name: worldcities
Geometry: None
Feature Count: 15493
Layer SRS WKT:
(unknown)
city: String (0.0)
city_ascii: String (0.0)
lat: String (0.0)
lng: String (0.0)
country: String (0.0)
iso2: String (0.0)
iso3: String (0.0)
admin_name: String (0.0)
capital: String (0.0)
population: String (0.0)
id: String (0.0)

ogr2ogr the Swiss army knife

Convert CSV to GeoPackage

In [10]:
%%time
!ogr2ogr -f GPKG worldcities.gpkg worldcities.csv \
  -oo X_POSSIBLE_NAMES=lng -oo Y_POSSIBLE_NAMES=lat -a_srs EPSG:4326
CPU times: user 9.59 ms, sys: 9.28 ms, total: 18.9 ms
Wall time: 887 ms

Use SQL statements

In [11]:
%%time
!ogr2ogr -f GPKG worldcities.gpkg worldcities.csv \
  -oo X_POSSIBLE_NAMES=lng -oo Y_POSSIBLE_NAMES=lat -a_srs EPSG:4326 \
  -sql "SELECT city, country, CAST(population AS integer) as population from worldcities where country = 'India'"
CPU times: user 4.59 ms, sys: 8.77 ms, total: 13.4 ms
Wall time: 381 ms

Rename layer in gpkg

In [12]:
%%time
!ogr2ogr -f GPKG worldcities.gpkg worldcities.csv \
  -oo X_POSSIBLE_NAMES=lng -oo Y_POSSIBLE_NAMES=lat -a_srs EPSG:4326 \
  -sql "SELECT city, country, CAST(population AS integer) as population from worldcities where country = 'India'" \
  -nln mycities
CPU times: user 4.89 ms, sys: 6.51 ms, total: 11.4 ms
Wall time: 443 ms
In [14]:
!ogrinfo -so -al worldcities.gpkg
INFO: Open of `worldcities.gpkg'
      using driver `GPKG' successful.

Layer name: mycities
Geometry: Point
Feature Count: 212
Extent: (69.670000, 8.180400) - (94.900000, 34.300000)
Layer SRS WKT:
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
FID Column = fid
Geometry Column = geom
city: String (0.0)
country: String (0.0)
population: Integer (0.0)

Exercise 6

In [16]:
!ogr2ogr mycities.shp worldcities.gpkg \
-t_srs EPSG:7755 -lco ENCODING=UTF-8 \
-lco SPATIAL_INDEX=Yes

Querying data using ogr2ogr

In [18]:
ls
1870_southern-india.jpg*           mycities.dbf
1870_southern-india.jpg.aux.xml    mycities.prj
5d16a93f1cf0f6000579ad2c.tif       mycities.qix
batch.py*                          mycities.shp
batch_parallel.py*                 mycities.shx
earth_at_night.jpg*                naip/
earth_at_night.tif                 precipitation.gpkg
earthquakes/                       prism/
gdal_stuff.qgz                     reservoirs.xml
geonames/                          reservoirs_india.tif
india-reprojected-polynomial.tif   reservoirs_india_proj.tif
india-reprojected-tps.tif          spatial_query.gpkg
india-reprojected-tps.tif.aux.xml  srtm/
india-with-gcp.vrt                 worldcities.csv*
landsat8/                          worldcities.gpkg
london_1m_dsm/                     worldcities.gpkg-shm
mycities.cpg                       worldcities.gpkg-wal
In [17]:
# Get country pop from city pop
!ogr2ogr country_pop.csv worldcities.gpkg \
  -sql "SELECT country, sum(population) as total_population from worldcities GROUP BY country"
ERROR 1: In ExecuteSQL(): sqlite3_prepare_v2(SELECT country, sum(population) as total_population from worldcities GROUP BY country):
  no such table: worldcities

Merging layers

In [19]:
cd earthquakes/
/Users/abharathi/Documents/gis_data/gdal-tools/earthquakes
In [20]:
ls
2020_01.geojson* 2020_04.geojson* 2020_07.geojson* 2020_10.geojson*
2020_02.geojson* 2020_05.geojson* 2020_08.geojson* 2020_11.geojson*
2020_03.geojson* 2020_06.geojson* 2020_09.geojson* 2020_12.geojson*
In [25]:
!which ogrmerge.py
/Users/abharathi/micromamba/envs/opengeo/bin/ogrmerge.py
In [26]:
%run -i /Users/abharathi/micromamba/envs/opengeo/bin/ogrmerge.py --help
ERROR: Unrecognized argument : --help
Usage: ogrmerge.py -o out_dsname src_dsname [src_dsname]*
            [-f format] [-single] [-nln layer_name_template]
            [-update | -overwrite_ds] [-append | -overwrite_layer]
            [-src_geom_type geom_type_name[,geom_type_name]*]
            [-dsco NAME=VALUE]* [-lco NAME=VALUE]*
            [-s_srs srs_def] [-t_srs srs_def | -a_srs srs_def]
            [-progress] [-skipfailures] [--help-general]

Options specific to -single:
            [-field_strategy FirstLayer|Union|Intersection]
            [-src_layer_field_name name]
            [-src_layer_field_content layer_name_template]

* layer_name_template can contain the following substituable variables:
     {AUTO_NAME}  : {DS_BASENAME}_{LAYER_NAME} if they are different
                    or {LAYER_NAME} if they are identical
     {DS_NAME}    : name of the source dataset
     {DS_BASENAME}: base name of the source dataset
     {DS_INDEX}   : index of the source dataset
     {LAYER_NAME} : name of the source layer
     {LAYER_INDEX}: index of the source layer
An exception has occurred, use %tb to see the full traceback.

SystemExit: 2
In [28]:
%run -i /Users/abharathi/micromamba/envs/opengeo/bin/ogrmerge.py -o earthquakes.gpkg *.geojson \
 -single -nln all_earthquakes -overwrite_ds

Write a new layer into the existing geopackage

In [32]:
!ogr2ogr earthquakes.gpkg earthquakes.gpkg \
-where "mag>4.5" -nln large_earthquakes -update

Spatial SQL

OGR supports dialect called SQLite dialect that is implemented by SpatiaLite. This is not as powerful as postGIS, but still useful.

ETL section

In [33]:
cd ..
/Users/abharathi/Documents/gis_data/gdal-tools
In [37]:
# %load batch.py
import os

input_dir = 'naip'

command = 'gdal_translate -of GTiff -co COMPRESS=JPEG {input} {output}'
for file in os.listdir(input_dir):
  if file.endswith('.jp2'):
    input = os.path.join(input_dir, file)
    filename = os.path.splitext(os.path.basename(file))[0]
    output =  os.path.join(input_dir, filename + '.tif')
    os.system(command.format(input=input, output=output))
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 5000, 5000
0...10...20...30...40...50...60...70...80...90...100 - done.
In [ ]:
# %load batch_parallel.py
import os
from multiprocessing import Pool
from timeit import default_timer as timer

input_dir = 'naip'

command = 'gdal_translate -of GTiff -co COMPRESS=JPEG {input} {output}'

def process(file):
    input = os.path.join(input_dir, file)
    filename = os.path.splitext(os.path.basename(file))[0]
    output =  os.path.join(input_dir, filename + '.tif')
    os.system(command.format(input=input, output=output))
    
files = [file for file in os.listdir(input_dir) if file.endswith('.jp2')]

if __name__ == '__main__':
  start = timer()
  p = Pool(4)
  p.map(process, files)
  end = timer()
  print(end - start)
  
  # start = timer()
  # for file in files:
  #   process(file)
  # end = timer()
  # print(end - start)
In [ ]: