How to parse parcel data

[momolomo-jpg]

I want to create a script that takes a KMZ line from google maps and creates a rectangle encompassing the line and uses this rectangle to select the parcel lines to store in a separate file. Thus far I have created a dwt file made from the dwg file that contains all the parcel lines in Washington state. I have also created a script that creates the rectangle in EPSG 4326 format. Now I'm trying to use geopandas to import the dwt and use bounding box and .intersect when reading the file to parse the lines. However, I couldn't get it to store any objects so I started troubleshooting by showing where the box was landing and it was way off its mark. I converted the dwt to EPSG 4326 but it won't work unless I import it in as 3857 first but even after changing to 4326 following this operation it places the geometry at 0,0 in matplot while the rectangle is way off. On top of that the time it takes between runs is fairly long. I'm not sure if the time it takes can be resolved but I'm not sure why the rectangle won't align when it shows up right in the KMZ.

[brendan-ward]

Not sure I follow what you are trying to do here.

Are the parcel data only available as DWG format, or are they also available in a more GIS-friendly format that definitely includes the correct CRS of parcel data? I haven't ever worked with DWG files, but it looks like maybe they don't have good support for CRS when read using GDAL: docs, which is what you are ultimately using via read_file in GeoPandas. To get around that, you'll need to set the CRS after reading the parcel data; do you know for certain that 3857 is the correct CRS of those source data? (it would be an odd choice for authoritative data like parcels)

Among other things, I'd suggest translating the DWG file into a more GIS friendly format (e.g., Flatgeobuff / GPKG are good options) via ogr2ogr and assign the correct CRS of the parcel data at that step. Then display the data in something like QGIS, or read out a tiny amount of records to display in the plot like you are doing in Python, to verify if they are displaying correctly.

In terms of slowness to run, it could be from a few things:

• reading the DWG as a format may be slow, esp. for data of this volume (no idea, never used DWG)
• using Fiona rather than pyogrio engine in read_file will be slower; you didn't state what you are using here
• doing the intersection will be really really slow if you don't use a spatial index; you didn't state how you were doing this, but what you described sounds like an intersects check against ALL parcels, which is a massive number.

Make sure to check out sjoin for doing an index-accelerated spatial join (if you want this as a spatial join) or use the spatial index query operation if you are trying to manage that more directly. Either way, you need to build a spatial index of the parcels, and then query that spatial index using the bounding rectangle of your lines from KML.

However, since you are dealing with a massive parcel dataset (>3M features), reading those data into memory and then building a spatial index will take some time and quite a lot of memory.

An alternative way to go about that would be to 1) convert the data into a GIS format that has a built-in spatial index (Flatgeobuff, GPKG, etc), and 2) select out only the data you need using the bbox parameter as part of read_file, which should then use the underlying on-disk spatial index to help find the possibly-intersecting features. You should then do a second pass over those data after reading to extract out exactly those ones that intersect, since using bbox may give you a few false positives.

[momolomo-jpg]

So the data is a DWG that I export to DXF. The correct CRS is 4326 but when I set it as 4326 I get an error (ValueError: aspect must be finite and positive) when using plot so I brought it in as 3857 then convert it into 4326 and that works for some reason. I'll try and find a good way to convert the DXF/DWG to a more GIS friendly format and check back in. The only issue I have with this is that I need to put it back into AutoCad DWG format so I was hoping to avoid another projection problem later.

[brendan-ward]

FYI - It looks like DWG is not supported for writing in GDAL, but DXF is; see: GDAL vector drivers.

[momolomo-jpg]

I figured out the problem was that I had the wrong EPSG I found a shape file and a geojson and when I imported them the region was different so I changed it in the DXF version and it worked. I still think that using the GIS version will be important though because it should save time which will be crucial when I place the parcels back in. (I was working on a reduced version for times sake)