1

I know that this is not an uncommon question, but I've struggled to find an answer that works for my application.

I'm trying to find the number of nearby permits for any given property. In order to do this, I've used a SpatiaLite index to create a buffer around each point and check the other points for intersection with that buffer.

The problem is that for some reason--regardless of the CRS that the parent layer of points is set to--the buffer will not project in meters or feet (depending on the projection). As you can see, I've set the CRS when creating the layer (and checked that it was correct after layer creation).

filepath = "filepath/permits.csv"
uri = filepath+"?delimiter=%s&xField=%s&yField=%s&crs=%s" % (",", "Latitude", "Longitude", "EPSG:2264") 
vlayer = QgsVectorLayer(uri, "Permits", "delimitedtext")


But the buffer still seems to be using decimal degrees for its projection. Any tips on how to get the buffer layer to project in feet or meters?

import sys

from qgis.core import *
from qgis.gui import *
from qgis.utils import iface

from PyQt4.QtGui import QApplication
from PyQt4.QtCore import QVariant

import itertools

# Initialize Application
app = QApplication([] , False)
QgsApplication.setPrefixPath("/Applications/QGIS.app/Contents/MacOS", True)
QgsApplication.initQgis()

# Create vector point layer from CSV
filepath = "filepath/permits.csv"
uri = filepath+"?delimiter=%s&xField=%s&yField=%s&crs=%s" % (",", "Latitude", "Longitude", "EPSG:3857") 
vlayer = QgsVectorLayer(uri, "Permits", "delimitedtext")
print vlayer.crs().authid()

#Use loaded vlayer as active vlayer to set the CRS
QgsMapLayerRegistry.instance().addMapLayer(vlayer)

# ID of each feature
_NAME_FIELD = 'PermitNum'

# Names of the new field to be added
_NEW_NEIGHBORS_FIELD = 'WithinMile'

# Create 2 new fields in the vlayer that will hold the # of neighbors
vlayer.startEditing()
vlayer.dataProvider().addAttributes([QgsField(_NEW_NEIGHBORS_FIELD, QVariant.String)])
vlayer.updateFields()

# Create a dictionary of all features
feature_dict = {f.id(): f for f in vlayer.getFeatures()}
print vlayer.crs().authid() 

# Build a spatial index
index = QgsSpatialIndex()
for f in feature_dict.values():
    index.insertFeature(f)

# Loop through all features and find features that touch each feature
for number, f in enumerate(feature_dict.values()):
    # Check to see if geometry is NoneType
    if f.geometry() is not None:
        print 'Working on %s' %  f[_NAME_FIELD]
        # Create buffer of given distance
        geom = f.geometry().buffer(5280,-1)
        intersecting_ids = index.intersects(geom.boundingBox())
        neighbors = len(intersecting_ids)
        print neighbors

vlayer.commitChanges()

print 'Processing complete.' 
  • Is the data truly in EPSG:2264? Or are the values in decimal degrees and labeled as being in NAD83 North Carolina state plane zone? – mkennedy Feb 26 '18 at 23:47
  • @mkennedy So it's a problem with the imported CSV? If that's the case, what's the simplest way to fully convert to the specified CRS? – stdmn Feb 26 '18 at 23:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.