trolleway_gis_snippets.md

Dockerfiles

Standart dockerfile

FROM registry.nextgis.com/sshd:0.1.0
#там ubuntu с gdal

ARG DEBIAN_FRONTEND=noninteractive
ARG APT_KEY_DONT_WARN_ON_DANGEROUS_USAGE=DontWarn

RUN apt-get update && apt-get install --no-install-recommends -y mc git nano wget tree
ARG uid=1000
ARG gid=1000
RUN groupadd -g $gid trolleway && useradd --home /home/trolleway -u $uid -g $gid trolleway  \
  && mkdir -p /home/trolleway && chown -R trolleway:trolleway /home/trolleway
RUN echo 'trolleway:user' | chpasswd

#у меня в деревне такой инет, что сразу все зависимости не выкачиваются, и этот уровень завершается.
#попробую ставить зависимости по частям, чтоб меньше качать
RUN apt-get install --no-install-recommends -y proj-data
RUN apt-get install --no-install-recommends -y python3-numpy
RUN apt-get install --no-install-recommends -y gdal-bin

RUN apt-get install --no-install-recommends -y python3-pip
#RUN apt-get install --no-install-recommends -y python3-psycopg2
RUN apt-get install --no-install-recommends -y time
RUN apt-get install --no-install-recommends -y osmctools osmium-tool aria2
RUN pip3 install tqdm
RUN pip3 install pyyaml
RUN pip3 install moz_sql_parser

#for ods2qlm

RUN apt update -y

#add to sudoers
RUN apt-get install -y apt-utils
RUN apt-get install -y sudo
RUN adduser trolleway sudo
RUN usermod -aG sudo trolleway

WORKDIR /data
#-p "$(openssl passwd -1 trolleway)"

Python

Remove directory with files

        directory='tmp'
        if os.path.exists(directory):
            import shutil
            shutil.rmtree('/'+tmp)

Create directory if not exists

        directory='tmp'
        if not os.path.exists(directory):
            os.makedirs(directory)

File walking in python

import os

for dirpath, dnames, fnames in os.walk("./"):
    for f in fnames:
        if f.endswith(".x"):
            x(os.path.join(dirpath, f))
        elif f.endswith(".xc"):
            xc(os.path.join(dirpath,f))

Connect to PostGIS

class Processor:

    conn_string = "host=192.168.250.1 user=trolleway dbname=osm_ch3 password="
     
   
    # get a connection, if a connect cannot be made an exception will be raised here
    conn = psycopg2.connect(conn_string)
    conn.autocommit = True #для vaccuum 
    # conn.cursor will return a cursor object, you can use this cursor to perform queries
    cursor = conn.cursor()
    

Import GeoJSON to PostGIS

    def import_local_geodata(self):
        cmd='ogr2ogr -f PostgreSQL PG:"{conn_string}" blacklist.geojson -nln blacklist -nlt MultiPolygon -overwrite'.format(conn_string=self.conn_string)
        os.system(cmd)

Import CSV to Python dict variables

    def read_csv(self,filename):
        print 'import tags table'
        print filename

        data=[]

        import csv
        with open(filename, 'rb') as f:
            reader = csv.reader(f)
            tags_list = list(reader)

        linecnt=0
        for line in tags_list:
            linecnt+=1
            if ((line[0] != '') and (linecnt>1)):
                element={}
                element['objectid']=line[0]
                element['name_ru']=line[4]
                element['sentinel_entity_id']=line[23]
                element['sentinel_earthexplorer_metadata_page']=line[23]
                element['landsat-8_product_identifier']=line[23]
                element['landsat-8_scene_identifier']=line[23]
                data.append(element)

        self.cities = data
        print self.cities

read_csv('Города и сцены - ndvi cities.csv')

Import CSV to PostgreSQL

 def import_book2(self):
        filename = '../../docs/railway/Консультант ТР4 Книга 2/book2_list.csv'

        #import csv to postgis

        print 'import Книга 2'
        sql = '''
                DROP TABLE IF EXISTS railgraph_book2 CASCADE; 
        CREATE TABLE railgraph_book2 (NAME varchar,CODES varchar,ROAD varchar,PAGE varchar,ROW varchar,TRANSIT varchar,CODE varchar);
        '''
        self.cursor.execute(sql)
        self.conn.commit()
        
        import csv
        with open(filename, 'rb') as f:
            reader = csv.reader(f)
            tags_list = list(reader)

        linecnt=0
        for line in tags_list:
            linecnt+=1
            if ((line[0] != '') and (linecnt>1)):
                print line[0],line[1]
                sql = '''INSERT INTO railgraph_book2 (NAME,CODES,ROAD,PAGE,ROW,TRANSIT,CODE) VALUES (
                '{col0}',
                '{col1}',
                '{col2}',
                '{col3}',
                '{col4}',
                '{col5}',
                '{col6}'
                )'''.format(
                    col0=line[0].replace("'","''"),
                    col1=line[1],
                    col2=line[2],
                    col3=line[3],
                    col4=line[4],
                    col5=line[5],
                    col6=line[6],
                    )
                print sql
                self.cursor.execute(sql)
                self.conn.commit()

PostgreSQL query in Python as dict

host=config.host
dbname=config.dbname
user=config.user
password=config.password

#Define our connection string
conn_string = "host={host} dbname={dbname} user={user} password={password}".format(host=host,dbname=dbname,user=user,password=password)
psql_string = '--host {host} -U {user} -d {dbname}'.format(host=host,dbname=dbname,user=user,password=password)

os.system('export PGPASS='+password)
 
conn = psycopg2.connect(conn_string)
conn.autocommit = True #для vaccuum
cursor = conn.cursor(cursor_factory=psycopg2.extras.DictCursor)

#Чтение таблицы из БД

sql = 'SELECT fid, address FROM suppliers'
cursor.execute(sql)
conn.commit()
rows = cursor.fetchall()
for row in rows:
    print dict(row)
    print '-'

Export PostGIS to CSV

    def generate_csv(self, region_code,filename):

        sql='''
SELECT 

ST_X(ST_Centroid(t1.wkb_geometry)) AS longitude,
ST_Y(ST_Centroid(t1.wkb_geometry)) AS latitude

'''
        sql=sql.format(region_code=region_code)


        self.cursor.execute(sql)
        self.conn.commit()
        rows = self.cursor.fetchall()
        records=[]            
        for row in rows:
            records.append(row)

        if os.path.exists(filename):
            os.remove(filename)

        with open(filename, 'w') as f:
            f.write('''"name_rus";"name_eng";"addr:district";"country_id"'''+"\n")

        with open(filename, 'a') as f:
            writer = csv.writer(f, delimiter=';',escapechar='"',quoting=csv.QUOTE_ALL)
           
            for row in records:
                writer.writerow(row)
        
        

Simpler export Postgresql to CSV using psycopg and COPY statement

       sqlv2 = '''
SELECT * from
        '''
       
        outputquery = "COPY ({0}) TO STDOUT WITH CSV".format(sqlv2)

        with open('_psycopg_result.csv', 'w') as f:
                self.cursor.copy_expert(outputquery, f)
		
	#add custom complex header to csv
	cat csv_header2.txt 170_v5_l1.csv> 170_v5_l2.csv	

Export PostGIS to geojson

def postgis2geojson(table):
    if os.path.exists(table+'.geojson'):
        os.remove(table+'.geojson')

        cmd='ogr2ogr -f GeoJSON spaces.geojson PG:"{conn_string}" -sql "SELECT name, area_full, area_r FROM spaces ORDER BY name" -nln areas -nlt MultiPolygon -overwrite'.format(conn_string=self.conn_string)
        os.system(cmd)
    print cmd
    os.system(cmd)

Escape SQL string in python

.replace("'","''")

Parse string

i = geofilename.rindex('_')
j = geofilename.rindex('.')
num = geofilename[i+2: j]

BeautifulSoup Parse table from web page

    def get_sentinel_vendor_id(self,url='https://earthexplorer.usgs.gov/metadata/10880/1731217/'):
        '''
        parse earthexplorer metadata page, return Vendor Product ID parameter
        '''

        s = requests.session()
        req = s.get(url)
        soup = BeautifulSoup(req.text, 'html.parser')
        data = []
        table = soup.find('table', attrs={'id':'metadataTable'})
        table_body = table.find('tbody')

        rows = table_body.find_all('tr')
        for row in rows:
            cols = row.find_all('td')
            cols = [ele.text.strip() for ele in cols]
            data.append([ele for ele in cols if ele]) # Get rid of empty values


        for row in data:
            if row[0] == 'Vendor Product ID':
                vendor_product_id=row[1]
        
        return vendor_product_id

PostGIS

LIKE operator

An underscore (_) in pattern stands for (matches) any single character; a percent sign (%) matches any string of zero or more characters.

Join attributes by location

SELECT
* 
FROM 
point p1, polygon p2
WHERE 
ST_Within(p1.geom, p2.geom)

Convert string to integer and return 0 if no digits

CAST((COALESCE(NULLIF(REGEXP_REPLACE(dist, '[^0-9]+', '', 'g'), ''), '0')) AS INTEGER)

Add primary key

ALTER TABLE test1 ADD COLUMN id SERIAL PRIMARY KEY;

Выбрать дорожную сеть как 1 обьект

SELECT 
     ST_Multi(ST_Collect(f.wkb_geometry)) as wkb_geometry
   FROM (SELECT (ST_Dump(wkb_geometry)).geom As wkb_geometry
        FROM
        highway_line WHERE highway='trunk') As f

Convert polygon layer to lines

DROP TABLE IF EXISTS ttemp;

CREATE TABLE ttemp AS
SELECT 
     ST_Boundary(ST_Multi(ST_Collect(f.wkb_geometry))) as wkb_geometry
   FROM (SELECT (ST_Dump(wkb_geometry)).geom As wkb_geometry
        FROM
        sea_basins ) As f

Select intersections of linestring and polygon layer

DROP TABLE IF EXISTS split_nodes;

CREATE TABLE split_nodes AS(
WITH cut_lines AS (SELECT 
     ST_Boundary(ST_Multi(ST_Collect(sub_area.wkb_geometry))) as wkb_geometry
   FROM (SELECT (ST_Dump(wkb_geometry)).geom As wkb_geometry
        FROM
        sea_basins) as sub_area)
SELECT
ST_Intersection(t1.wkb_geometry, t2.wkb_geometry) as wkb_geometry
FROM sea_routes_split as t1, cut_lines AS t2
WHERE ST_INTERSECTS(t1.wkb_geometry, t2.wkb_geometry)
GROUP BY ST_Intersection(t1.wkb_geometry, t2.wkb_geometry)
)

In case of error "Operations on mixed SRID"

UPDATE isaa_points SET wkb_geometry=ST_SetSRID(wkb_geometry,4326)

Snap point layer to line layer

UPDATE client_points
SET    is_snapped=true,
       wkb_geometry=ST_ClosestPoint(network.wkb_geometry, client_points.wkb_geometry)
FROM   (
              SELECT ST_Multi(St_collect(f.wkb_geometry)) AS wkb_geometry
              FROM   (
                            SELECT (ST_dump(wkb_geometry)).geom as wkb_geometry
                            FROM   highway_line
                            WHERE  highway IN ('trunk', 'motorway' ) AS f ) AS network
WHERE    st_distance(st_closestpoint(network.wkb_geometry, client_points.wkb_geometry)::geography,client_points.wkb_geometry::geography) < 30000

Generate points at roads intersections

--Generate points at roads intersections
DROP TABLE IF EXISTS lines_intersections ;
CREATE TABLE lines_intersections AS
SELECT      
    ST_Intersection(a.wkb_geometry, b.wkb_geometry) AS wkb_geometry,
    Count(Distinct a.ogc_fid)
FROM
    sea_routes as a,
    sea_routes as b
WHERE
    ST_Touches(a.wkb_geometry, b.wkb_geometry) OR ST_Crosses(a.wkb_geometry, b.wkb_geometry) 
    --touches returns T-shape crossings, crosses returns X-shape crossings
    AND a.ogc_fid != b.ogc_fid
GROUP BY
    ST_Intersection(a.wkb_geometry, b.wkb_geometry)
;

Generate fake service roads from points to nearest highway

--table of points on roads near to terminal. Points generating at all hoghways in radius.
    -- таблица точек на дорогах напротив терминалов. Генерируются точки на всех ближайших дорогах в заданом радиусе
CREATE TEMP TABLE cut_points AS (
SELECT 
terminals.id AS id, line_id AS line_id,
ST_ClosestPoint(subquery.wkb_geometry::geography, terminals.geom::geography) AS wkb_geometry,
ST_Distance(terminals.geom::geography, ST_ClosestPoint(subquery.wkb_geometry, terminals.geom)::geography) AS dist
FROM    (SELECT 
            ST_Multi(St_collect(network.wkb_geometry)) AS wkb_geometry,
            network.line_id
        FROM   
            (SELECT (ST_dump(wkb_geometry)).geom as wkb_geometry, id AS line_id FROM roads  ) AS network
        GROUP BY network.line_id    ) AS subquery,
        terminals 
WHERE    
    st_distance(st_closestpoint(subquery.wkb_geometry, terminals.geom)::geography,terminals.geom::geography) < $MAX_CONN_LEN
);
COMMENT ON TABLE cut_points IS 'nodes at network near poinst. Here points at all roads in distance';

-- Table for select only one nearest point for each terminal
 -- таблица для выборки только одной ближайшей точки от каждого терминала до дороги
CREATE TEMP TABLE nearest_cutpoints AS
  (SELECT cut_points.id,
          cut_points.line_id,
          cut_points.wkb_geometry
   FROM cut_points,

     (SELECT id,
             min(dist) AS dist
      FROM cut_points
      GROUP BY id) AS subquery
   WHERE cut_points.id=subquery.id
     AND cut_points.dist=subquery.dist);
COMMENT ON TABLE nearest_cutpoints IS 'Select only one nearest network node for each point';

-- Connections of terminals to roads
    -- таблица связей терминалов с дорогами
  DROP TABLE IF EXISTS conn;
  CREATE TABLE conn AS (
    SELECT 
    terminals.id AS id,
    ST_MakeLine(terminals.geom,nearest_cutpoints.wkb_geometry) AS wkb_geometry,
    ST_Distance(terminals.geom::geography, nearest_cutpoints.wkb_geometry::geography) AS dist
    FROM terminals JOIN nearest_cutpoints ON terminals.id = nearest_cutpoints.id
  );
COMMENT ON TABLE conn IS 'Lines from points to network';

--Insert nodes to network at intercections with lines from points. Snap road to point
--Добавление точку в дороги в местах примыкания пепендикуляров от поставщиков
UPDATE roads
SET wkb_geometry=ST_Snap(roads.wkb_geometry,points_near_way.wkb_geometry,$TOLERANCE)
FROM
  (SELECT ST_Collect(points.wkb_geometry) AS wkb_geometry,
          points.line_id
   FROM nearest_cutpoints AS points
   GROUP BY points.line_id
    ) AS points_near_way
WHERE points_near_way.line_id=roads.id;

planet_osm_rels to simple table

DROP TABLE IF EXISTS unnesting_with_rn;
CREATE TEMPORARY TABLE unnesting_with_rn AS 
	select id, unnest, row_number() OVER (PARTITION BY id)  FROM 
	(
	SELECT 
	id, 
	unnest(members) 
	FROM
	planet_osm_rels
	) AS unnesting;

SELECT members.id as relation_id, members.unnest as member,  roles.unnest as role, members.row_number,roles.row_number  FROM 
unnesting_with_rn AS members JOIN unnesting_with_rn AS roles ON roles.id=members.id and roles.row_number=members.row_number+1 
WHERE roles.row_number % 2 = 0;	 

#TODO: calc in 3857, not 4326

Суммарная протяжённость маршрутов автобуса, троллейбуса, трамвая в базе PostGIS импортированной из OSM

CREATE OR REPLACE FUNCTION unnest_rel_members_ways(ANYARRAY) RETURNS
SETOF ANYELEMENT
LANGUAGE SQL AS $$SELECT substring($1[i] from E'w(\\d+)') FROM
generate_series(array_lower($1,1),array_upper($1,1)) i WHERE 
$1[i] LIKE 'w%' /*only ways*/
AND /*exclude platforms*/
($1[i+1] ='' 
OR $1[i+1] IN ('forward','backward','highway')
)
;$$;



SELECT sum(ST_Length_Spheroid(ST_Transform(planet_osm_line.way,4326),'SPHEROID["WGS 84",6378137,298.257223563]'))/1000 
FROM
planet_osm_line,
(
SELECT

  unnest_rel_members_ways(members) as members 
from planet_osm_rels
WHERE 

tags::VARCHAR LIKE '%route,trolleybus%'
            OR tags::VARCHAR LIKE '%route,tram%'
            OR tags::VARCHAR LIKE '%route,bus%'
            OR tags::VARCHAR LIKE '%route,share_taxi%'
) AS unnestWays
WHERE unnestWays.members::bigint = planet_osm_line.osm_id

Location of shops near line

--Запрос выдаёт пикетные отметки по одиночной линии из таблицы trolleybus, на которых находятся магазины (из таблицы prodykty). Работает с EPSG:32637 (для Москвы)
SELаECT 
row_number() OVER (ORDER BY 1) AS i, produkty."@id",
produkty.geom,
round((ST_line_locate_point(trolleybus.geom,produkty.geom) * ST_Length(trolleybus.geom)/1000)::numeric ,1) AS pk,
COALESCE(produkty.name,produkty.operator,produkty."@id") AS name
FROM produkty, trolleybus
WHERE ST_Intersects(produkty.geom ,  ST_Buffer(trolleybus.geom,500))
ORDER BY pk

Select ways from OSM for bus line

DROP TABLE IF EXISTS route_lines ;
CREATE TEMPORARY TABLE route_lines AS 
SELECT *, members::hstore AS members_hstore,
tags::hstore->'ref' AS ref, 
tags::hstore->'operator' AS operator FROM planet_osm_rels
WHERE tags::hstore -> 'route' = 'trolleybus' AND tags::hstore->'ref' = '2';

SELECT * FROM route_lines;

DROP TABLE IF EXISTS target_route;
CREATE TABLE target_route AS
SELECT 
planet_osm_line.name,
planet_osm_line.way
FROM 
planet_osm_line JOIN route_lines ON (route_lines.members_hstore ? CONCAT('w',planet_osm_line.osm_id::varchar)  )
;
SELECT * FROM target_route;

Extract bus line from OSM PostGIS database to geojson

message() {
#Display text message in terminal
    dialog --backtitle 'Data processing' --title 'gnu dialog'  --nook --nocancel --infobox "$1" 4 40
    sleep 1
}

rm -rf temp
mkdir temp
message 'Calculate route'

cat > temp/query.sql <<- EOF

DROP TABLE IF EXISTS route_lines ;
CREATE TEMPORARY TABLE route_lines AS 
SELECT *, members::hstore AS members_hstore,
tags::hstore->'ref' AS ref, 
tags::hstore->'operator' AS operator FROM planet_osm_rels
WHERE tags::hstore -> 'route' = 'trolleybus' AND tags::hstore->'ref' = '2';

SELECT * FROM route_lines;

DROP TABLE IF EXISTS target_route;
CREATE TABLE target_route AS
SELECT 
planet_osm_line.name,
planet_osm_line.way
FROM 
planet_osm_line JOIN route_lines ON (route_lines.members_hstore ? CONCAT('w',planet_osm_line.osm_id::varchar)  )
;

EOF

psql  -d $dbname -a -f temp/query.sql
#sleep 1

rm -r target_line.geojson 
ogr2ogr  -f "GeoJSON"  "target_line.geojson" PG:"dbname=gis" target_route

message 'Data extracted to geojson'
#

Select dublicates

select code_from, code_to, count(*)
from railgraph_book1_edges
group by code_from, code_to
HAVING count(*) > 1

Select rows not in another table

SELECT name 
FROM   railgraph_edges l 
WHERE  NOT EXISTS (
   SELECT i.code              -- it's mostly irrelevant what you put here
   FROM   railway_stations i
   WHERE  l.code_to = i.code
   );

Buffer with dissolve results

CREATE TEMPORARY TABLE buffers100 ON COMMIT DROP AS 
SELECT st_buffer(ST_Transform(wkb_geometry,32637),100) AS wkb_geometry,100 AS buffer_size FROM moscow_crossings WHERE ways_cnt = 4;

WITH
 clusters(wkb_geometry) AS 
      (SELECT ST_CollectionExtract(unnest(ST_ClusterIntersecting(wkb_geometry)), 3) 
         FROM buffers100),
 multis(id, wkb_geometry) AS 
      (SELECT row_number() over() as id, wkb_geometry FROM clusters)
 SELECT ST_UNION(wkb_geometry) FROM 
      (SELECT id, (ST_DUMP(wkb_geometry)).geom AS wkb_geometry FROM multis) d GROUP BY id;

Delete dublicates

DELETE FROM osm2esr
WHERE ogc_fid IN (SELECT ogc_fid
              FROM (SELECT ogc_fid,name,
                             ROW_NUMBER() OVER (partition BY esr ORDER BY ogc_fid) AS rnum
                     FROM osm2esr) t
              WHERE t.rnum > 1 ORDER BY name);
	      

2 approach: http://technobytz.com/most-useful-postgresql-commands.html

Update with join

UPDATE vehicles_vehicle AS v 
SET price = s.price_per_vehicle
FROM shipments_shipment AS s
WHERE v.shipment_id = s.id 

Select lines crosses with state boundary

       DROP TABLE IF EXISTS test ;
        CREATE  TABLE test AS
SELECT railgraph_book1_edges.id ,railgraph_book1_edges.way
FROM railgraph_book1_edges, state_boundary_railways 
WHERE ST_Intersects(railgraph_book1_edges.way, ST_Boundary(state_boundary_railways.wkb_geometry));

Split line layer at features at intersections

Withouth attributes May fail if one line crosses other many times https://gis.stackexchange.com/questions/71270/split-line-at-intersection-and-attach-attributes

DROP TABLE IF EXISTS sea_routes_split;

CREATE TABLE sea_routes_split AS(
select  (st_dump(st_union(wkb_geometry))).geom AS wkb_geometry 
from sea_routes);
ALTER TABLE sea_routes_split ADD COLUMN id SERIAL PRIMARY KEY;

Print CREATE TABLE from psql to display withouth comments

pg_dump --dbname=postgresql://trolleway:trolleway@127.0.0.1:5432/gis \
--schema-only \
--format plain \
--clean \
--no-privileges \
--table suppliers\
| sed -e '/^--/d' \
| cat -s

Import OSM PBF to PostGIS

        sql='''
        DROP TABLE IF EXISTS planet_osm_line CASCADE;
        DROP TABLE IF EXISTS planet_osm_point CASCADE;
        DROP TABLE IF EXISTS planet_osm_polygon CASCADE;
        DROP TABLE IF EXISTS planet_osm_roads CASCADE;

        '''
        self.cursor.execute(sql)
        self.conn.commit()

        os.system('export PGPASS='+self.pgpass)
    
        print 'pbf to o5m'
        cmd='osmconvert {filename}.osm.pbf -o={filename}.o5m'.format(filename=filename)
        print cmd        
        os.system(cmd)

        print 'o5m tag filtration'
        cmd='osmfilter {filename}.o5m --drop-author --keep="{fl}"   --out-o5m >{filename}-filtered.o5m'.format(filename=filename, fl=self.generate_filter_string())
        print cmd        
        os.system(cmd)

        print 'o5m to pbf'
        cmd='osmconvert {filename}-filtered.o5m -o={filename}-filtered.pbf'.format(filename=filename)
        print cmd        
        os.system(cmd)


        print 'pbf to postgis'
        cmd='osm2pgsql {osm2pgsql_config} --password --slim  --create --multi-geometry --latlon   -C 12000 --number-processes 3 --style special.style {filename}-filtered.pbf'.
format(osm2pgsql_config=config.osm2pgsql,filename=filename)
        print cmd        
        os.system(cmd)

Extract buildings from OSM dump

export PGPASS=
filename=france-latest
filter="building=*"
osmconvert $filename.osm.pbf -o=$filename.o5m
osmfilter $filename.o5m --drop-author --keep="$filter"   --out-o5m >$filename-filtered.o5m
osmconvert $filename-filtered.o5m -o=$filename-filtered.pbf

Openstreetmap

Download and clip pbf

CLIP="RU-CHU.parts.geojson"
DUMP="far-eastern-fed-district-latest.osm.pbf"

# deploy
wget --output-document=ogr2poly.py https://trac.openstreetmap.org/browser/subversion/applications/utils/osm-extract/polygons/ogr2poly.py?format=txt

#run
wget --timestamping https://download.geofabrik.de/russia/far-eastern-fed-district-latest.osm.pbf

mv $CLIP clip.geojson
python ogr2poly.py clip.geojson

#ogr2poly produce one file for each feature
#osm.pbf extension is required for JOSM (https://wiki.openstreetmap.org/wiki/JOSM/Plugins/PBF)
osmconvert $DUMP -B=clip_0.poly -o=clipped_0.o5m
osmconvert $DUMP -B=clip_1.poly -o=clipped_1.o5m

DATETIME=$(date +%Y-%m-%d_%H-%M-%S)
osmconvert clipped_0.o5m clipped_1.o5m -o=$DATETIME.osm.pbf

Filter foot routing graph from OSM

url="https://download.geofabrik.de/russia/crimean-fed-district-latest.osm.pbf"
src="crimean-fed-district-latest.osm.pbf"
wget --timestamping $url

osmconvert $src -b=33.9814,44.9218,34.2314,45.052 -o=clipped.osm
#foot graph 
osmfilter clipped.osm --keep="highway=" --drop="access=private access=no foot=no highway=motorway highway=construction highway=proposed highway=path highway=track highway=cycleway highway=raceway" >streets.osm 

Import OSM to PostGIS

#assume you has database gis and allow rin to postgis withouth password

filename='.osm.pbf'
osm2pgsql --create --slim --latlon  --database gis $filename

Import OSM to PostGIS with relation_members table

#https://medium.com/@trolleway/import-osm-pbf-file-to-postgis-database-with-relation-members-25b84fc10208

#Install latest osmosis
#current version in Ubuntu fails while read pbf
wget http://bretth.dev.openstreetmap.org/osmosis-build/osmosis-latest.tgz
mkdir osmosis
mv osmosis-latest.tgz osmosis
cd osmosis
tar xvfz osmosis-latest.tgz
rm osmosis-latest.tgz
#End of osmosis install
psql -d gis -c "CREATE extension IF NOT EXISTS hstore"
psql -d gis -f /home/trolleway/osmosis/script/pgsnapshot_schema_0.6.sql
osmosis/bin/osmosis \

--read-pbf file="filtered.pbf" \
--write-pgsql host="localhost" database="gis" user="" password=""

Ubuntu

Zip all files in current directiry

zip -r extract.zip .

Generate cendroids in ogr2ogr

time ogr2ogr -progress -overwrite -lco ENCODING=UTF-8  centroid.shp ../building-polygon.shp -dialect sqlite -sql "SELECT ST_Centroid(geometry), building,  name from \"building-polygon\" ORDER BY A_HSNMBR LIMIT 1000000"

Merge all csv into Shapefile

#!/bin/bash
consolidated_file="./union.shp"
for i in $(find . -name '*.csv'); do
    if [ ! -f "$consolidated_file" ]; then
        # first file - create the consolidated output file
        ogr2ogr -f "ESRI Shapefile" -lco encoding="UTF-8" -oo X_POSSIBLE_NAMES=Lon* -oo Y_POSSIBLE_NAMES=Lat* $consolidated_file $i
    else
        # update the output file with new file content
        ogr2ogr -f "ESRI Shapefile" -update -append  -lco encoding="UTF-8" -oo X_POSSIBLE_NAMES=Lon* -oo Y_POSSIBLE_NAMES=Lat* $consolidated_file $i
    fi
done

Filter pbf by tags

filename=central-fed-district-latest
osmconvert $filename.osm.pbf -o=$filename.o5m
osmfilter $filename.o5m --drop-author --keep="all highway="   --out-o5m >$filename-filtered.o5m
osmconvert $filename-filtered.o5m -o=$filename-filtered.pbf

Start Opera Use Opera VPN

overpass-turbo.eu

/*
This has been generated by the overpass-turbo wizard.
The original search was:
“admin_level=4 and name = "Орловская область"”
*/
[out:json][timeout:25];
// gather results
(
  // query part for: “admin_level=4 and name="Орловская область"”

  relation["admin_level"="4"]["name"="Орловская область"]({{bbox}});
);
// print results
out body;
>;
out skel qt;

Open geojson in NextGIS QGIS

export to poly using plugin
rename poly to border.poly

sourcefile=central-fed-district-latest
wget http://download.geofabrik.de/russia/$sourcefile.osm.pbf
osmconvert $sourcefile.osm.pbf -o=temp.o5m
osmconvert temp.o5m -B=border.poly --complete-multipolygons -o=result.osm

Merge geodata and convert to pbf for routing

# Assume you have layers with same crosing points
# frist git clone ogr2osm
# see http://wiki.openstreetmap.org/wiki/Ogr2osm

ogr2ogr -progress -overwrite -geomfield wkb_geometry -t_srs EPSG:4326 -f "ESRI Shapefile"  merged_lines.shp 30.geojson  
ogr2ogr -progress -update  -append -geomfield wkb_geometry -t_srs EPSG:4326 -f "ESRI Shapefile" -nlt LINESTRING -where "OGR_GEOMETRY='LineString'"  merged_lines.shp 36.geojson 


#добавление атрибутов
ogrinfo merged_lines.shp -sql "ALTER TABLE merged_lines ADD COLUMN highway character(25)"
ogrinfo merged_lines.shp -dialect SQLite -sql "UPDATE merged_lines SET highway = 'primary'"


#Конвертация из ogr в osm
python /home/trolleway/GIS/GIS/project112_asgu/ogr2osm/ogr2osm.py --force --id 999999999 --positive-id merged_lines.shp 

#Конвертация osm в pbf для компактности
osmconvert merged_lines.osm -o=merged_lines.osm.pbf

Resize photos (bash)

mkdir resized
time for f in `find . -name "*.jpg"` ;do     convert $f -resize 50% -quality 55% -sampling-factor 4:4:4 resized/$f.resized.jpg; done

Batch reproject vector layers

Reproject to EPSG:4326

for i in *.shp do; ogr2ogr -t_srs EPSG:4326 reprojected/$i $i; done 

Reproject to МСК-00

rm -r reprojected 
mkdir reprojected

for i in *.shp; do ogr2ogr -t_srs "+proj=tmerc +lat_0=0 +lon_0=46.05 +k=1 +x_0=1300000 +y_0=-4714743.504 +ellps=krass +towgs84=23.57,-140.95,-79.8,0,0.35,0.79,-0.22 +units=m +no_defs" reprojected/$i $i; done 

Get GeoJSON for OSM feature by id

pip install osm2geojson

wget -O result.osm "https://overpass-api.de/api/interpreter?data=[out:xml][timeout:25];(relation(2697545););out meta;>;out meta qt;" && osm2geojson -f  result.osm result.geojson

Windows batch video

Extract frames from video and make timelapse (Windows terminal)

MKDIR %~p1frames
DEL /Q %~p1frames
ffmpeg -i "%1" -r 2 -f image2 -vcodec mjpeg -qscale 1 "%~p1frames/output_%%05d.jpg"
ffmpeg -f image2 -i "%~p1frames/output_%%05d.jpg" -qscale 1 -b 3000k -y %~p1%~n1_speed%~x1

Extract frames from video (Windows terminal)

MKDIR frames
DEL /Q frames
ffmpeg -i ".MOV" -r 2 -f image2 -vcodec mjpeg -qscale 1 "frames/output_%05d.jpg"

More quality extract frames from video for photogrammetry (Windows terminal)

MKDIR frames
DEL /Q frames
ffmpeg -i ".mp4" -r 2 -f image2 -vf fps=2 -qscale 1 "frames/output_%05d.png"

Make a gif from frames (Unix)

ffmpeg -f image2  -pattern_type glob -i '*.png' -framerate 1 -qscale 1 animation.gif

Convert MOV to MP4 Windows

ffmpeg -i movie.mov -vcodec copy -acodec copy out.mp4

Convert MOV to MP4 and drop audio Windows

ffmpeg -i movie.mov -vcodec copy -an out.mp4

Create video from frames

ffmpeg  -framerate 60    -f image2 -i "Untitled%04d.jpg" output.mp4

Cut video in windows

ffmpeg -ss 00:04:07 -i kiev.mp4   -t 17 -c copy kievcrop.mp4

Extract frames from some scenes of film as jpg

Paste in Windows console

SET video = "Tusjacha.okon (1967).avi"
SET folder = "Tyshacha_Okon_1967_Minsk1"
MKDIR %folder%
ffmpeg -ss 00:03:00 -i %video% -t 12 -r 2 -f image2 -vcodec mjpeg -qscale 1 "%folder%/frame_%05d.jpg"

ffmpeg trim mp3 audio

ffmpeg -i "Jean-Michel Jarre - Calypso.mp3" -acodec copy  -ss 00:01:58 -t 20 -y   music.mp3

ffmpeg video with static image and music

ffmpeg -loop 1 -y -i ueit.png -i sound.mp3 -acodec copy   -vcodec libx264 -t 1622 "Bob Cobert - Music From Express To Terror (1979).avi"

ffmpeg convert video from left chanel sound only to stereo

ffmpeg -i "e:\Android\data\com.cyberlink.powerdirector.DRA140225_01\files\paznis_Full HD.mp4" -vcodec copy -af "pan=mono|c0=c0" pazniy_mono.mp4

ffmpeg join video

(echo file 'e:\el_photos\Xiaomi_360\20180515\rotated\VID_20180515_224000AA_5528113_veka.MP4' & echo file 'e:\el_photos\Xiaomi_360\20180515\rotated\VID_20180515_224233AA_5528113_veka.MP4' & echo file 'e:\el_photos\Xiaomi_360\20180515\rotated\VID_20180515_224427AA_5528113_veka.MP4' )>list.txt
ffmpeg -safe 0 -f concat -i list.txt -c copy -vcodec copy -vsync drop -r 60 join_veka.mp4

git clone https://github.com/google/spatial-media.git
python spatial-media\spatialmedia  -i join_veka.mp4 join_veka_injected.mp4
del join_veka.mp4

Agisoft Photoscan

Change photo path in Agisoft Photoscan

{projectname}.files\0\0\frame.zip\doc.xml
<frame>
  <cameras>
    <camera camera_id="12">
      <photo path="