You can cast your geometry as XY before outputting it to GeoJSON:
ogr2ogr -f csv -dialect sqlite -t_srs epsg:3857 -sql "select
AsGeoJSON(CastToXY(geometry)) AS geom, * from PVS_19_v1_vtd_21031"
partnership_shapefiles_19v1_21031/PVS_19_v1_vtd_21031.shp -dim 2
One approach is to use the SQLite SQL dialect. It seems to work at least with a point shapefile where I digitized a few points, some of them in the same location.
ogrinfo -dialect sqlite -sql "select geometry, count(*) from duplicate_points group by geometry" duplicate_points.shp
Layer name: SELECT
Geometry: Unknown (any)
Feature Count: 4
Complete workaround for converting each feature to an individual layer and then rasterize each one it can be observed in following code (where my own paths to layers were used):
from osgeo import gdal
from osgeo import ogr
from osgeo import gdalconst
raster_file = "/home/zeito/pyqgis_data/mask.tif"
shape_file = "/home/zeito/pyqgis_data/test2.shp"
If you have proj 6 installed, you can query the proj.db file to get the list of EPSG codes it supports. On my machine, the following works:
sqlite3 /usr/share/proj/proj.db "SELECT code FROM projected_crs WHERE auth_name = 'EPSG';"
Your code is perfect.
I think you are not working in a projection in meters e.g mollweide.
All you need to do is reproject the layer to -- +proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m no_defs
Then run the code again and it would be 100%