The two functions you're comparing don't produce the same output. gdal_polygonize will create polygons connecting equal-valued pixels (in a single polygon). Compared to the QGIS function, which converts every pixel to a unique point/polygon regardless it's value.
It's unclear to me why your result of gdal_polygonize appears to contain points instead of polygons. That shouldn't be possible.
You could try doing the same from Python with gdal.Polygonize. When I use the snippet below, it seems to work as I would expect. The output are the red-lines in the image below, overlayed on the input raster.
from osgeo import gdal, ogr, osr
ras_file = 'NT60ne_DTM_2m.tif'
vec_file = 'NT60ne_DTM_2m.gpkg'
ds_ras = gdal.OpenEx(ras_file)
src_band = ds_ras.GetRasterBand(1)
srs_wkt = ds_ras.GetProjection()
srs = osr.SpatialReference()
srs.ImportFromWkt(srs_wkt)
drv = ogr.GetDriverByName('GPKG')
ds_vec = drv.CreateDataSource(vec_file)
lyr = ds_vec.CreateLayer('', srs=srs, geom_type=ogr.wkbPolygon)
lyr.CreateField(ogr.FieldDefn("elevation", ogr.OFTReal))
fld = lyr.GetLayerDefn().GetFieldIndex("elevation")
gdal.Polygonize(src_band, None, lyr, fld)
srs = None
ds_ras = None
ds_vec = None
A subset, showing only the lower-right corner:
edit:
The snippet below is an example of writing each pixel to a feature (Point/Polygon) to disk, with the raster value attached as an attribute. You could avoid using Shapely & Geopandas if needed, but it would require a bit more code to use only OGR to create the vector-file.
The same technique could be use to output a Polygon but the creation of the WKT would be slightly more complicated, requiring the extent of each pixel. That can be created by offsetting the currently used center coordinates with +/- half the resolution.
A method like this, using string formatting/parsing for a large number of points, will never be extremely performant.
import geopandas as gpd
from shapely import wkt
import pandas as pd
import numpy as np
from osgeo import gdal, ogr, osr
ras_file = 'NT60ne_DTM_2m.tif'
vec_file = 'NT60ne_DTM_2m_points.gpkg'
# read raster data & metadata
ds = gdal.OpenEx(ras_file)
gt = ds.GetGeoTransform()
xs = ds.RasterXSize
ys = ds.RasterYSize
srs_wkt = ds.GetProjection()
nodata = ds.GetRasterBand(1).GetNoDataValue()
data = ds.ReadAsArray()
ds = None
# mask to ignore nodata
assert np.isfinite(nodata) # else use np.isfinite / np.isnan
valid_data = data != nodata
# (outer) extent of raster
ulx, xres, _, uly, _, yres = gt
lrx = ulx + xs * xres
lry = uly + ys * yres
# convert extent to center of pixel
ulx = ulx + xres/2
lrx = lrx - xres/2
uly = uly + yres/2
lry = lry - yres/2
# create coordinate grids
mapy, mapx = np.mgrid[uly:lry+yres:yres, ulx:lrx+xres:xres]
# create normal DataFrame
df = pd.DataFrame(dict(mapx=mapx[valid_data], mapy=mapy[valid_data], elevation=data[valid_data]))
# Use WKT Point to add the geometry (just a string column)
df["geometry"] = df.apply(lambda row: wkt.loads(f"POINT ({row.mapx} {row.mapy})"), axis=1)
# convert to actual GeoDataFrame
df = gpd.GeoDataFrame(df, geometry='geometry', crs=srs_wkt)
# save to disk
df.to_file(vec_file, driver="GPKG")
