Writing compressed vector data for access in QGIS

Question

I'm trying to find an optimum vector file format to be used to effectively store data (compressed binary?) and select subregions.

I have a quite large ascii point file containing 3 columns (x and y projected coordinates and the point ID). A test file is about 20MB, containing approximately 600000 points, but I plan to create even larger files with significantly more points. After writing those files to ascii, I'm searching for a way (i.e. file format) to store the data more effectively while still allowing for the following two tasks:

• the files shall be loaded into QGIS (v2.14.17)
• subregions shall be selected (e.g. via GDAL/OGR v2.1 "ogr2ogr -spat") and rewritten in the original ascii format

Do you have any suggestions, which file format shall be used here?

EDIT: Just to add some more information the data looks like

CoordX,CoordY,ID  
1200400,-399600,1
1201200,-399600,2
1202000,-399600,3
1202800,-399600,4
...

I've created a VRT file as follows:

<OGRVRTDataSource>
    <OGRVRTLayer name="points">
        <SrcDataSource>points.csv</SrcDataSource>
        <GeometryType>wkbPoint</GeometryType>
        <LayerSRS>EPSG:3031</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="CoordX" y="CoordY"/>
        <Field name="X" src="CoordX" type="integer" />
        <Field name="Y" src="CoordY" type="integer" />
        <Field name="ID" src="ID" type="integer" />
    </OGRVRTLayer>
</OGRVRTDataSource>

As Luke suggested (What format to use for highly compressed vector data?) I tried SpatiaLite with compressed geometry: ogr2ogr -f SQLite -dsco SPATIALITE=YES -lco COMPRESS_GEOM=YES -lco COMPRESS_COLUMNS=x,y,id points.sqlite points.vrt but the result was a 90MB file (compared to 14MB of the now slightly modified CSV-Input).

Playing around a bit showed, that leaving away the -lco COMPRESS_* did not have any effect on the 90MB file size. Instead, leaving away -dsco SPATIALITE=YES resulted in a 26MB file (I know, this is not SpatiaLite any more). I'm working on Ubuntu 14.04. GDAL 2.1.0 was installed using the UbuntuGIS-PPA. Is there any additional package (besides libspatialite5 v4.1.1) necessary to enable compression in SpatiaLite?

Answer 1

Great AndreJ, that did the trick! I was working with NetCDF for a long time but never thought about using it for vector data.
ogr2ogr -f NetCDF points.nc points.shp produced a NetCDF file which was still 17MB, i.e. more than the input. But this NetCDF was indeed possible to be loaded into QGIS. I didn't find any option to tell ogr2ogr to use compression within the netcdf. However, now that I know the structure how a NetCDF should look to be recognized as vector data by GDAL/OGR, I wrote a Python script to read the data and rewrite it as NetCDF similar to the ogr2ogr-output but including compression. It worked like a charm!! Now the the ~600000 points just occupy 111kB. It needs some time to load in QGIS (I guess it is converted back into a temporary shapefile in the background) but finally it works!

#!/usr/bin/python
import sys
import os
import re
from netCDF4 import Dataset

### read argument 
if(len(sys.argv) < 2 or not re.search("\.csv$",sys.argv[1])):
    print "Usage: ",sys.argv[0],"<input>.csv"
    exit(1)

input_csv=sys.argv[1]
basnam=os.path.splitext(os.path.basename(input_csv))[0]
output_nc="%s.nc"%(basnam)

### read data file
(x,y,id)=([],[],[])
fcsv = open(input_csv,'r')
lines = fcsv.readlines()
for line in lines:
    cols = line.split(',')
    if(cols[0]=="CoordX"): continue
    x.append(int(cols[0]))
    y.append(int(cols[1]))
    id.append(int(cols[2]))
fcsv.close()
print "read",len(id),"records"

### Writing NetCDF file
nc_out = Dataset(output_nc, "w", format="NETCDF4_CLASSIC")
nc_out.Conventions = "CF-1.6"
nc_out.ogr_layer_name = basnam
nc_out.featureType = "point"
## dimension
nc_out.createDimension('record', None)

## variables
# coordinates
nc_out_var = nc_out.createVariable('x','f8',('record'),zlib=True)
nc_out_var.standard_name = "projection_x_coordinate"
nc_out_var.units = "m"
nc_out.variables['x'][:]    = x
nc_out_var = nc_out.createVariable('y','f8',('record'),zlib=True)
nc_out_var.standard_name = "projection_y_coordinate"
nc_out_var.units = "m"
nc_out.variables['y'][:]    = y

# projection
nc_out_var = nc_out.createVariable('polar_stereographic','c')
nc_out_var.grid_mapping_name = "polar_stereographic" ;
nc_out_var.straight_vertical_longitude_from_pole = 0. ;
nc_out_var.false_easting = 0. ;
nc_out_var.false_northing = 0. ;
nc_out_var.latitude_of_projection_origin = -90. ;
nc_out_var.standard_parallel = -71. ;
nc_out_var.long_name = "CRS definition" ;
nc_out_var.longitude_of_prime_meridian = 0. ;
nc_out_var.semi_major_axis = 6378137. ;
nc_out_var.inverse_flattening = 298.257223563 ;
nc_out_var.spatial_ref = "PROJCS[\"WGS_84_Antarctic_Polar_Stereographic\",GEOGCS[\"GCS_WGS_1984\",DATUM[\"WGS_1984\",SPHEROID[\"WGS_84\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.017453292519943295]],PROJECTION[\"Polar_Stereographic\"],PARAMETER[\"latitude_of_origin\",-71],PARAMETER[\"central_meridian\",0],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1]]" ;

# fields
nc_out_var = nc_out.createVariable('CoordX','i4',('record'),zlib=True)
nc_out_var.long_name = "Field X"
nc_out_var.grid_mapping = "polar_stereographic"
nc_out.variables['CoordX'][:]    = x

nc_out_var = nc_out.createVariable('CoordY','i4',('record'),zlib=True)
nc_out_var.long_name = "Field Y"
nc_out_var.grid_mapping = "polar_stereographic"
nc_out.variables['CoordY'][:]    = y

nc_out_var = nc_out.createVariable('ID','i4',('record'),zlib=True)
nc_out_var.long_name = "Field ID"
nc_out_var.grid_mapping = "polar_stereographic"
nc_out.variables['ID'][:]    = id

nc_out.close()

print "Finished writing",output_nc

I consider the topic as solved, however an option to create a compressed NetCDF using ogr2ogr would be desirable.