How to Convert Between Geospatial Formats: Shapefile, GeoPackage, GeoJSON and Cloud Optimized GeoTIFF

Data arrives in whatever format its maker preferred. You then convert it to the one that suits your tool, your colleague or your website. Converting is simple. The hard part is knowing what each format cannot keep.

This guide uses the Ludhiana district boundary and the clipped DEM from the guides on clipping and reprojecting. The data zip is ludhiana_workflow_data.zip. Unzip it next to your notebook.

The vector formats

For your own work, GeoPackage is the best default. Use Shapefile only when someone else's tool needs it, and GeoJSON when the data is going to a web map.

Vector, route A: QGIS

1. Add the layer with Layer > Add Layer > Add Vector Layer.
2. Right-click it in the Layers panel and choose Export > Save Features As.
3. Choose the Format: GeoPackage, ESRI Shapefile or GeoJSON.
4. Choose the File name. For a GeoPackage, give the Layer name as well, because one file can hold many layers.
5. Check the CRS. For GeoJSON, choose EPSG:4326, as explained below.
6. Under Select fields to export, untick the fields you do not need.
7. Click OK. QGIS adds the new layer to the project.

The Processing tool Save vector features to file does the same without the dialog, and it can run for a batch of layers.

Vector, route B: Python

You need Python with geopandas, rasterio and numpy.

import geopandas as gpd

# Step 1: read the district boundary. GeoPandas picks the right driver from the file.
district = gpd.read_file("ludhiana_workflow_data/ludhiana_district.geojson")
print("fields:", district.columns.tolist())
print("CRS: EPSG", district.crs.to_epsg(), "| features:", len(district))

# Step 2: add two fields with long names, to see what each format does with them
district["district_name_long"] = "Ludhiana"
district["population_estimate"] = 1000

fields: ['shapeName', 'shapeISO', 'shapeID', 'shapeGroup', 'shapeType', 'geometry']
CRS: EPSG 4326 | features: 1

Now write the same layer in all three formats.

import os

# Step 1: GeoPackage. One file. The layer name goes inside it.
district.to_file("district.gpkg", layer="district", driver="GPKG")

# Step 2: Shapefile. Several files with the same name.
district.to_file("district.shp")

# Step 3: GeoJSON. One text file.
district.to_file("district.geojson", driver="GeoJSON")

# Step 4: what was written, and how big is it
for name in sorted(os.listdir(".")):
    if name.startswith("district."):
        print(f"{name:18s} {os.path.getsize(name):7d} bytes")

district.cpg             5 bytes
district.dbf           757 bytes
district.geojson      9298 bytes
district.gpkg       102400 bytes
district.prj           145 bytes
district.shp          6124 bytes
district.shx           108 bytes

The Shapefile is five files, and they must stay together: .shp holds the shapes, .dbf the attributes, .shx an index, .prj the CRS and .cpg the text encoding. If you send only the .shp, the receiver gets shapes with no attributes and no CRS.

What each format kept

# Step 1: read each file back and compare the field names and the CRS
for path in ["district.gpkg", "district.shp", "district.geojson"]:
    back = gpd.read_file(path)
    print(path, "| EPSG", back.crs.to_epsg(), "| features:", len(back))
    print("   fields:", [c for c in back.columns if c != "geometry"])

district.gpkg | EPSG 4326 | features: 1
   fields: ['shapeName', 'shapeISO', 'shapeID', 'shapeGroup', 'shapeType', 'district_name_long', 'population_estimate']
district.shp | EPSG 4326 | features: 1
   fields: ['shapeName', 'shapeISO', 'shapeID', 'shapeGroup', 'shapeType', 'district_n', 'population']
district.geojson | EPSG 4326 | features: 1
   fields: ['shapeName', 'shapeISO', 'shapeID', 'shapeGroup', 'shapeType', 'district_name_long', 'population_estimate']

Look at the field names in the Shapefile. district_name_long became district_n and population_estimate became population. A Shapefile limits a field name to 10 characters, and the writer cuts the extra ones without a warning that you will notice. If two long names start with the same ten letters, they collide. The GeoPackage and the GeoJSON kept both names whole. This is the most common reason to leave the Shapefile behind.

GeoJSON is always in longitude and latitude

The GeoJSON standard says coordinates are longitude and latitude in WGS 84 (EPSG:4326). A reader assumes that. Data in a UTM zone, written to GeoJSON, can be read in the wrong place by a tool that follows the standard. Reproject to EPSG:4326 before you write GeoJSON.

# Step 1: the district in a UTM zone (metres)
district_utm = district.to_crs(32643)
print("UTM coordinates of the first corner:", [round(v) for v in district_utm.total_bounds[:2]])

# Step 2: for GeoJSON, go back to longitude and latitude first
district_utm.to_crs(4326).to_file("district_web.geojson", driver="GeoJSON")
back = gpd.read_file("district_web.geojson")
print("written for the web: EPSG", back.crs.to_epsg(), "| first corner:", [round(float(v), 3) for v in back.total_bounds[:2]])

UTM coordinates of the first corner: [534023, 3381191]
written for the web: EPSG 4326 | first corner: [75.356, 30.56]

A GeoPackage or a Shapefile can keep any CRS, and it should be the CRS you analyse in. Convert to EPSG:4326 only for the web.

Check that nothing was lost

# Step 1: compare the original with the GeoPackage round trip. Area is measured in a projected CRS.
original_area = district.to_crs(32643).area.sum() / 1e6
roundtrip_area = gpd.read_file("district.gpkg").to_crs(32643).area.sum() / 1e6
print("area before:", round(original_area, 1), "km2 | after:", round(roundtrip_area, 1), "km2")
print("same number of features:", len(district) == len(gpd.read_file("district.gpkg")))

area before: 3696.5 km2 | after: 3696.5 km2
same number of features: True

After any conversion, check three things: the number of features, the CRS and the field names. For polygons, add the total area.

Raster: GeoTIFF and Cloud Optimized GeoTIFF

A GeoTIFF is the default raster format. A Cloud Optimized GeoTIFF, or COG, is a GeoTIFF whose data is arranged in a way that lets a program read a small window or a zoomed-out view without reading the rest of the file. It has two parts that make that possible: the data is stored in tiles, which are small square blocks, and the file holds overviews, which are smaller copies of the raster at half, quarter and eighth size. A COG is still an ordinary GeoTIFF that any tool can open. That is why the DEM and the Sentinel-2 files in the download guides could be read from a server one window at a time.

Make a COG when a raster will be put on a server, shared online, or read by a web map. For files that stay on your own disk, an ordinary compressed GeoTIFF is fine.

Raster, route A: QGIS

QGIS builds a COG in two steps. In the Translate dialog, the option for the output format cannot be changed, so the COG layout comes from creation options.

1. Add dem_ludhiana_clip.tif. Work on a copy, because the next step changes the file.
2. Open Raster > Miscellaneous > Build Overviews (Pyramids). Set Overview levels to 2 4 8 16. Set Overviews format to Internal (if possible). Click Run. The overviews are stored inside the file.
3. Open Raster > Conversion > Translate (Convert Format).
4. Set Input layer to the copy with overviews.
5. Under Additional creation options, add three options: COMPRESS=DEFLATE, TILED=YES and COPY_SRC_OVERVIEWS=YES.
6. Save the result as dem_ludhiana_cog.tif and click Run.

Do not type -of COG into the Additional command-line parameters. QGIS already adds the output format of its own, and the command fails with the message "Duplicate argument -of". The three creation options above give the same layout.

Open the properties of the new file. The Information page shows a block size of 256 by 256, which means it is tiled, and the overviews are listed.

Raster, route B: Python

import rasterio
import rasterio.shutil

# Step 1: the clipped DEM, as it is now. Is it tiled, and does it have overviews?
with rasterio.open("ludhiana_workflow_data/dem_ludhiana_clip.tif") as src:
    print("original: tiled:", src.profile.get("tiled"), "| block:", src.block_shapes[0], "| overviews:", src.overviews(1))

# Step 2: write a COG. The COG driver adds the tiles and the overviews itself.
rasterio.shutil.copy("ludhiana_workflow_data/dem_ludhiana_clip.tif", "dem_ludhiana_cog.tif", driver="COG", compress="deflate")

# Step 3: check the result
with rasterio.open("dem_ludhiana_cog.tif") as cog:
    print("cog: tiled:", cog.profile.get("tiled"), "| block:", cog.block_shapes[0], "| overviews:", cog.overviews(1), "| compression:", cog.compression.value)
    print("same size:", cog.width, "x", cog.height, "| CRS EPSG", cog.crs.to_epsg(), "| NoData", cog.nodata)

original: tiled: True | block: (256, 256) | overviews: []
cog: tiled: True | block: (512, 512) | overviews: [2, 4, 8] | compression: DEFLATE
same size: 3497 x 1665 | CRS EPSG 4326 | NoData -9999.0

The COG has overviews at 2, 4 and 8 times smaller, and its data is in tiles of 512 by 512 cells. The original had no overviews. The overviews are the reason a COG is a little larger than the plain file.

# Step 1: do the heights agree? Read both and compare.
with rasterio.open("ludhiana_workflow_data/dem_ludhiana_clip.tif") as a, rasterio.open("dem_ludhiana_cog.tif") as b:
    same = bool((a.read(1) == b.read(1)).all())
    print("identical values:", same)

identical values: True

A conversion to a COG changes the layout and not the values. The values must be identical.

Which format for what

Common mistakes

• Sending a Shapefile as a single .shp file.
• Not noticing that a Shapefile cut the field names.
• Writing GeoJSON in a UTM zone.
• Typing -of COG into the Translate dialog in QGIS.
• Building a COG from a raster that had no NoData value set.
• Not checking the feature count and the CRS after a conversion.

Key takeaways