python - Arcpy remove invalid character

I've built a model that takes the name of the Excel worksheet containing xy data, and uses that name as the name of the output shapefile. ( Worksheet "massachusetts_sightings" in an xlsx becomes "massachusetts_sightings$.shp"). I do this using Parse Path. Notice that a "$" is added.

If I then feed this shapefile, massachusetts_sightings$.shp into a clip tool (say, to clip out points in the data that are not in my study area), I get an error for an invalid character in the clip tool. My suspicion is the $. I know I could remove this with Calculate Value...but I can't figure out the code block.

My parse path output is called Value:

Can someone help troubleshoot? I haven't used Python in years-

# Import arcpy module
import arcpy


# Load required toolboxes
arcpy.ImportToolbox("Model Functions")

# Set Geoprocessing environments
arcpy.env.scratchWorkspace = "X:\\Tools\\XY\\point_upload\\Scratch"
arcpy.env.workspace = "X:\\Tools\\XY\\point_upload\\Data"

# Script arguments
XY_Table = arcpy.GetParameterAsText(0)


Define_Coordinate_System = arcpy.GetParameterAsText(1)

Projection = arcpy.GetParameterAsText(2)

Geographic_Transformation = arcpy.GetParameterAsText(3)

Study_Area = arcpy.GetParameterAsText(4)

Projected_and_Clipped_Data = arcpy.GetParameterAsText(5)
if Projected_and_Clipped_Data == '#' or not Projected_and_Clipped_Data:

    Projected_and_Clipped_Data = "X:\\Tools\\XY\\point_upload\\Data\\%fixed%_Clip.shp" # provide a default value if unspecified

# Local variables:
xy_event = XY_Table
v_defined = xy_event
defined_data = v_defined
v_projected = defined_data
Value = XY_Table

# Process: Parse Path

arcpy.ParsePath_mb(XY_Table, "NAME")

# Process: Make XY Event Layer
tempEnvironment0 = gp.scratchWorkspace
arcpy.env.scratchWorkspace = "X:\\Tools\\XY\\point_upload\\Scratch"
tempEnvironment1 = gp.workspace
arcpy.env.workspace = "X:\\Tools\\XY\\point_upload\\Data"
arcpy.MakeXYEventLayer_management(XY_Table, "", "", xy_event, "", "")
arcpy.env.scratchWorkspace = tempEnvironment0
arcpy.env.workspace = tempEnvironment1


# Process: Copy Features
tempEnvironment0 = gp.scratchWorkspace
arcpy.env.scratchWorkspace = "X:\\Tools\\XY\\point_upload\\Scratch"
arcpy.CopyFeatures_management(xy_event, v_defined, "", "0", "0", "0")
arcpy.env.scratchWorkspace = tempEnvironment0

# Process: Define Projection
arcpy.DefineProjection_management(v_defined, Define_Coordinate_System)


# Process: Project
tempEnvironment0 = gp.scratchWorkspace
arcpy.env.scratchWorkspace = "X:\\Tools\\XY\\point_upload\\Scratch"
tempEnvironment1 = gp.workspace
arcpy.env.workspace = "X:\\Tools\\XY\\point_upload\\Data"
arcpy.Project_management(defined_data, v_projected, Projection, Geographic_Transformation, "")
arcpy.env.scratchWorkspace = tempEnvironment0
arcpy.env.workspace = tempEnvironment1

# Process: Clip

arcpy.Clip_analysis(v_projected, Study_Area, Projected_and_Clipped_Data, "")

# Process: Calculate Value
arcpy.CalculateValue_management("\"%Value%\".strip(\"$\")", "\\n", "String")

Exporting Feature GeoJSON from PostGIS?

I'd like to do a GeoJSON dump of a PostGIS table, but I want to export features (the geometry and the properties) not just the geometry. I've been digging into the ST_AsGeoJSON function but it looks like I can only get the coordinate information (which makes sense since it's a geometry function).

For example:

select st_asgeojson(the_geom) from street_centerline limit 1;

Returns:

{"type":"MultiLineString","coordinates":[[[-65.591776562805038,41.682190576167052],[-65.591998971284028,41.682082119060382],[-65.592001213509064,41.682081025737766],[-65.593689871787177,41.681257533373952],[-65.595415661879244,41.680415888937219],[-65.595440519465640,41.680403765889309],[-65.595603134242481,41.680324459445771]]]}

I'm curious if anyone knows of a simple way to get the property information as well. Has anyone written a pgsql2geojson script yet?

qgis - How does Zonal Statistics work exactly?

I'm a bit confused on how the zonal statistics tools work in QGIS (either the 'Zonal Statistics' tool under the Raster tools from the QGIS geoalgorithms toolbox or the 'Raster statistics from polygons' tool under the SAGA Vector<->Raster tools).

Basically I have a raster map layer and some small polygons in a separate layer. I want to assign each polygon the value of the raster below. In case a polygon overlaps multiple raster cells, I just want to have the maximum.

Yet the output doesn't make any sense and I can't get behind the reasoning of those results (which are obviously not the maximum of the underlain raster cell):

Any ideas what I'm doing wrong or what I didn't consider yet? The results with the python API in qgis.analysis.QgsZonalStatistics() (according to this question here) are much better, but it does only offer count, mean and sum...but for some reason no max or min.

python - Clipping vs Intersecting

I was curious why clipping shapefiles was just so difficult and whether it was different to intersecting shapefiles?

You will see from one my questions that I gave up trying to clip a polyline with a polygon (here)

I tried various different methods:

1. 
   

   Conversion ogr2ogr (too slow)

   
   

   import subprocess
   subprocess.call(["C:\Program Files\GDAL\ogr2ogr", "-f", "ESRI Shapefile", "-clipsrc", clipping_shp, output_shp, input_shp])
   print('Done!')
   
   



2. 
   

   SQL ogr (too slow)

   
   

   from osgeo import ogr
   ogr.UseExceptions()
   ogr_ds = ogr.Open('K:/compdata/OS_Shape_Files/os_open_roads/trim', True)  # Windows: r'C:\path\to\data'
   start = time.clock()
   print('Starting:')
   print(ogr_ds)
   SQL = """\
   
   SELECT ST_Intersection(A.geometry, B.geometry) AS geometry, A.*, B.*
   FROM ROADLINK_Project A, cut_se59ef_poly_60 B
   WHERE ST_Intersects(A.geometry, B.geometry);"""
    layer = ogr_ds.ExecuteSQL(SQL, dialect='SQLITE')
   # copy result back to datasource as a new shapefile
   layer2 = ogr_ds.CopyLayer(layer, 'h1_buf_int_ct3')
   # save, close
   layer = layer2 = ogr_ds = None
   print("Finished in %.0f seconds" % time.clock() - start)
   




3. 
   

   pyclipper (couldn't get this to work)

   
   

   solution = pc.Execute2(pyclipper.CT_INTERSECTION,     pyclipper.PFT_NONZERO, pyclipper.PFT_NONZERO)
   



4. 
   

   rTree (couldn't get this work as I have python 3.4)

   
   

   import fiona
   from shapely.geometry import shape, mapping
   import rtree
   
   
   bufSHP = produce_poly_shape
   intSHP = 'K:/compdata/OS_Shape_Files/os_open_roads    /trim/ROADLINK_Project.shp'
   ctSHP  = 'output.shp'
   
   start = time.clock()
   print('Starting')
   
   with fiona.open(bufSHP, 'r') as layer1:
   with fiona.open(ctSHP, 'r') as layer2:
       # We copy schema and add the  new property for the new resulting shp
   
       schema = layer2.schema.copy()
       schema['properties']['uid'] = 'int:10'
       # We open a first empty shp to write new content from both others shp
       with fiona.open(intSHP, 'w', 'ESRI Shapefile', schema) as layer3:
           index = rtree.index.Index()
           for feat1 in layer1:
               fid = int(feat1['id'])
               geom1 = shape(feat1['geometry'])
               index.insert(fid, geom1.bounds)
   
   
           for feat2 in layer2:
               geom2 = shape(feat2['geometry'])
               for fid in list(index.intersection(geom2.bounds)):
                   if fid != int(feat2['id']):
                       feat1 = layer1[fid]
                       geom1 = shape(feat1['geometry'])
                       if geom1.intersects(geom2):
                           # We take attributes from ctSHP
                           props = feat2['properties']
                           # Then append the uid attribute we want from the other shp
   
                           props['uid'] = feat1['properties']['uid']
                           # Add the content to the right schema in the new shp
                           layer3.write({
                               'properties': props,
                               'geometry': mapping(geom1.intersection(geom2))
                           })
   
   print("Finished in %.0f seconds" % time.clock() - start)
   

Finally, I looked at this but couldn't get it to work with python 3.4.

Then in R I tried:

1. 
   

   gIntersection(x,y, byid=TRUE)

   
   


2. 
   

   gDifference

   
   


3. 
   
   

   custom gClip command I found online

   
   

However with a large polylines shapefile (around 500mb) and a very small polygon shapefile (1mb) the intersection would take a day. This led me to think perhaps I am doing something incredibly stupid and there is a quick clip command?

For example, in arcGIS the intersection command takes a few seconds so surely clipping is even easier? I know very little about GIS but to me it seems as simple as aligning the two layers by one coordinate (assuming same projection) and then selecting the outside of the clipper shape (e.g. in paint) and just deleting it from the other layer. I guess obviously this isn't the case ..

Thus my goal: I create a polygon in my code and I want to clip the UK roads network I load in to the extent of that polygon, buffer the points, dissolve and then output as an array containing the polygon coordinates.

I don't need to retain any features that an intersection would provide just a straight clip. I guess I also don't strictly need shape files and could convert my 500mb UK road network into a geodatabase?

qgis - Layer with custom UI not working properly

I have a problem with QGIS 2.4 when introducing layers with a custom UI. More precisely, widgets such as CheckBoxes are not recognised as such (appear as Text Edit in Attributes) and the corresponding Checked value appears as NULL (although checked by default in the UI file).

When I am manually changing the specific field:

instudy_id = layer.fieldNameIndex("instudy")
layer.setEditorWidgetV2(instudy_id, 'CheckBox') 

Everything is normal. The corresponding part of the UI file is:

Sorry if my explanation is poor, I am quite new to QGIS. Hope the screenshots help.

To conclude:

• 
  

  some options of my Widgets (defined in my .ui file) are overwritten by QGIS (e.g. "Checked" or "Unchecked" state of a QCheckBox)

  
  


• 
  

  others, which are introduced in the corresponding .py file (like QgsEditorWidgetWrapper.fromWidget( in study ).setEnabled(False) ) are not respected by QGIS.

  
  

What is the proper way to introduce a custom / customise an .ui file ? Is there anywhere detailed information/examples on the use of newer methods like QgsVectorLayer.setEditorWidgetV2, setEditorWidgetV2Config, QgsEditorWidgetWrapper ?

arcmap - Counting raster values in each polygon using ArcGIS Desktop?

I have a Watershed with ~100 smaller watersheds as a shapefile. I want to count the number of different land use values (they are in a 30x30m raster) in each of the small watersheds.

I tried combining them but nothing turns out they way I need it. Which tool or routine can I use to get this?

Answer

One of the powerful tools in ArcGIS is the Tabulate Area tool, which allows you to calculate the area of each land use within each polygon in the polygon shapefile (watershed in your case).

The output will be a table that computes the area of land use in each polygon. If you want to know the number of pixels of land use simply you can create a new column and divide the area by 900 (30 x 30) and you will get the number of pixels of land use within each polygon.

Spatial Analyst extension is required to use Tabulate Area tool.

style - Is Geoserver's "styler" plugin dead?

The "Styler" plugin for Geoserver promised flexible styling of your data, using the GeoExt library:

The user guide tells you to first install the REST plugin, and then extract the styler zip file into the /data/www Geoserver directory. However, the download page for Geoserver 2.2 does not include the REST plugin. Installing older versions (I'm using Geoserver 2.2) failed for me, with one error saying "This web site needs a different Google Maps API key" (from /geoserver/www/styler/ - I have a Google Maps API).

REST seems to have disappeared recently: versions up to 2.0.3 contained the option to download the REST zip package, whereas after that there is no mention of REST, and the installation guide remains.

I did find some evidence that REST may have died, however:

This feature is no longer updated as an independent GeoServer download, but is available as part of the OpenGeo Suite. For more information on its development see OpenGeo project page

If Geoserver's "styler" plugin really is dead, then I should probably change the question to something more general like "how do I make OpenGeo's Styler work on Geoserver?" Just wanted to check with others who've had similar experience.

Context I've set-up Geoserver on an AWS EC2 instance (testing map can be seen here). What I'm trying to do is to make those boring grey polygons colourful and interesting, based on their attributes. I've tried exporting the sld that discribes the style below into Geoserver, but, at present, this does not work.

Answer

The OpenGeo Suite's GeoExplorer now has the styling capability that Styler used to have. That must be why they shelved it.

You can get the code and installation instructions at Github.

qgis - Add raster value to features AttributeError: 'module' object has no attribute 'TYPE_RASTER'?

When I try to run the "add raster value to features" (available into SAGA Tools), I get this error:

Traceback (most recent call last):
  File "/usr/share/qgis/python/plugins/processing/gui/AlgorithmDialog.py", line 249, in accept
    msg = self.alg._checkParameterValuesBeforeExecuting()

  File "/usr/share/qgis/python/plugins/processing/core/GeoAlgorithm.py", line 239, in _checkParameterValuesBeforeExecuting
    return self.checkParameterValuesBeforeExecuting()
  File "/usr/share/qgis/python/plugins/processing/algs/saga/SagaAlgorithm.py", line 343, in checkParameterValuesBeforeExecuting
    param.datatype == dataobjects.TYPE_RASTER):
AttributeError: 'module' object has no attribute 'TYPE_RASTER'

I'm using QGIS 2.18.10 and Python 2.7.12 on Ubuntu 16.04

I don't have any clue to solve this error.

Getting shapefile of river from OpenStreetMap?

I do not want to download the whole map of an area. I just want rivers and other waterbodies present in an area in shapefiles so that I can use them in QGIS. I tried geofabrik but it is giving whole map and that also of full country while i need data of a city. Right now I am trying JOSM. will update you when it will works.

Answer

Use the Overpass Turbo API!

First, I'd go to openstreetmap.org,

You can either

1) Search for the feature you're after. I chose "south platte river", which runs through Denver. This gives the fields and tags that are used by OSM to store the data:

2) Identify the tags and values of the features you're after by

1. Zooming all the way into the map


2. Click on the layers icon on the right (the three sheets of paper)


3. Click on the last menu entry (Map data or something similar in your language)


4. The features on the map turn blue (make sure you're zoomed in far enough to see


5. Click on the feature you're after



6. The Tags and Values appear on left side of the screen, and you can proceed below...

Then go on over to the Overpass Turbo page, then click Wizard

Using information, the name value is South Platte River, and the waterway value is river, so you can build a query like this:

name="South Platte River" and waterway=river

Then click "build and run query"

The query will run and the result will show on the map:

Next click the "Export" option:

I like the geoJSON option

Click "Save"

Open the file in QGIS, and away you go!

You can do a 'save as' to save it as a new type of vector layer...

In your case, you could also use the waterway=river query to get all the rivers in the area you're after, and you can draw a manual selection box to narrow down the geography.

@underdark showed me this.

qgis - How to merge many individual KML layers into one?

I'd like to open a USGS KML file of bedrock geology for an entire county, but there are over 2,000 individual polygon layers. Does anyone know if there is a way to merge all of these layers together into one file before opening in QGIS?

coordinate system - Should point data be equidistant projected when using ArcGIS IDW spatial interpolation?

I am wondering if anyone can clarify whether point data should be projected to an equidistant projection when using ArcGIS IDW spatial interpolation?

I am working on a dataset from Western North America spanning about 30 degrees of latitude. The data are currently in Lat/Long (NAD83).

Does ArcGIS "project on the fly" or somehow adjust for latitude when calculating distances to my sample points during the interpolation procedure or should I be supplying everything in a projection that preserves distances?

Answer

IDW works by finding the data points located nearest each point of interpolation, weighting the data values according to a given power p of the distances to those points, and forming the weighted average. (Often p = -2.)

Suppose there is some amount of distance distortion around an interpolation point that is the same in all directions. This will multiply all distances by some constant value x. The weights therefore all get multiplied by x^p. Because this does not change the relative weights, the weighted average is the same as before.

When the distance distortion changes with direction, this invariance no longer holds: data points in some directions now appear (on the map) relatively closer than they should while other points appear relatively further. This changes the weights and therefore affects the IDW predictions.

Consequently, for IDW interpolation we would want to use a projection that creates roughly equal distortions in all directions from each point on the map. Such a projection is known as conformal. Conformal projections include those based on the Mercator (including Transverse Mercator (TM)), Lambert Conic, and even Stereographic.

It is important to realize that conformality is a "local" property. This means that the distance distortion is constant across all bearings only within small neighborhoods of each point. For larger neighborhoods involving greater distances, all bets are off (in general). A common--and extreme--example is the Mercator projection, which is conformal everywhere (except at the poles, where it is not defined). Its distance distortion becomes infinite at sufficiently large north-south distances from the Equator, while along the Equator itself it's perfectly accurate.

The amount of distortion in some projections can change so rapidly from point to point that even conformality will not save us when the nearest neighbors are far from each other or near the extremes of the projection's domain. It is wise, then, to choose a conformal projection adapted to the study region: this means the study region is included within an area where its distortion is the smallest. Examples include the Mercator near the Equator, TM along north-south lines, and Stereographic near either pole. In the conterminous US, the Lambert Conformal Conic is often a good default choice when the reference latitudes are placed within the study region but near its northern and southern extremes.

These considerations usually are important only for study regions that extend across large countries or more. Within small countries or states of the US, popular conventional coordinate systems exist (such as various national grids and State Plane coordinates) which introduce little distance distortion within those particular countries or states. They are good default choices for most analytical work.

python - Fiona - Preffered method for defining a schema

I'm fairly new to the world of python and GIS - so this is definitely a beginner question.

I am following Tom Macwrights small guide to Shapely and Fiona, and I noticed that even when opening the previously created shapefile, it doesn't contain the schema I defined using:

schema = { 
      'geometry': 'Point', 
      'properties': { 'name': 'str' } }

This is evident because when I try to call it using schema = input.schema.copy(), it spits out AttributeError: 'function' object has no attribute 'schema'

Reading the Fiona documentation, it states that:

the schema of its record type (a vector file has a single type of record, remember) is accessed via a read-only schema attribute.

Does that mean that it is not possible to store a schema in the shapefile,a nd therefore I'm always required set up the variables in the beginning, rather than reading from the file using input.schema.copy()? Or do I setup the schema using Records?

If so, how do you read from a csv file (as in Tom's example) and place the data into the records to be read from input.schema.copy()? is it helpful to do this, or am I just complicating things and should define the schema every time i open the shapefile?

Thanks for the help - I'm really interested in developing my understanding of the processes occurring.

Answer

I don't really understand your problem. Some explanations: for example, we want to

1. create a new shapefile


2. 
   

   modify the original schema: for that, it's easier to just copy things to a new shapefile and make the changes as you copy

   
   
   
   
   • Create a new Polyline shapefile:
   
   
   
   
   

import fiona
# schema: it is a simple dictionary with geometry and properties as keys
schema = {'geometry': 'LineString','properties': {'test': 'int'}}
# for defining the geometry, you need Shapely
from shapely.geometry import LineString, mapping
# two simples geometries
lines = [LineString([(272830.63,155125.73),(273770.32,155467.75)]),LineString([(273536.47,155914.07),(272033.12,152265.71)])]

with fiona.open('myshp.shp', 'w', 'ESRI Shapefile', schema) as layer:
    for line in lines:
        # filling schema
        elem = {}
        # geometry with mapping function of shapely
        elem['geometry'] = mapping(line) 
        # attribute value (the same here)
        elem['properties'] = {'test': 145}
        # writing element in the file
        layer.write(elem)

• Now we want to modify the schema of the original shapefile in a new shapefile:

1) Open the original shapefile:

shapefile =fiona.open('myshp.shp')
#read the schema
schema2 = shapefile.schema
print schema2
{'geometry': 'LineString', 'properties': OrderedDict([(u'test', 'int:10')])}

2) As it is a dictionary, it is easy to add new fields/keys in the properties:

schema2['properties']['string']='str'

3) Now we create a new shapefile copying myshp.shp with the new schema :

with fiona.open('myshp.shp', 'r') as input:
    schema = schema2
    # writing the new shapefile
    with fiona.open('myshp_copy.shp', 'w', 'ESRI Shapefile', schema) as output:
        for elem in input:

            # add the new attribute value
            elem['properties']['string']="hello"
            output.write({'properties': elem['properties'],'geometry': mapping(shape(elem['geometry']))})

• verification

c = fiona.open('myshp_copy.shp')
c.schema
{'geometry': 'LineString', 'properties': OrderedDict([(u'test', 'int:10'), (u'string', 'str')])}

# first element of the shapefile
c.next()
{'geometry': {'type': 'LineString', 'coordinates': [(272830.63, 155125.73000000001), (273770.32000000001, 155467.75)]}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'test', 145), (u'string', u'hello')])}

• Conclusion

If you don't want to modify the shapefile, it is easier with schema.copy() that is only used to get a copy of the original schema (no definition here)

with fiona.open('myshp.shp', 'r') as input:
    schema = input.schema.copy()

     with fiona.open('myshp_copy2.shp', 'w', 'ESRI Shapefile', schema) as output:
         for elem in input:
             output.write({'properties': elem['properties'],'geometry': mapping(shape(elem['geometry']))})

arcmap - Creating relationship database in ArcGIS Desktop

I am creating a file geodatabase with a feature class. I have created attribute columns in which a few of them have domain value YES/NO but when I select, NO, I have to provide a measurement value as per the condition.

How do I design such a database?

The objective of this database is to collect field information through the ArcGIS collector App where the field workers will measure the length of the roadside drainage if it does not satisfy the length specification.

QGIS Python script for creating project file

Is there a way to create QGIS project file(.qgs) using only Python script outside QGIS Python console. I am trying to publish a map using QGIS Web Client using python.

For this, I am trying to do below

1. Create .qgs file dynamically


2. Add Legend and Layers to the .qgs file and save


3. Set WFS,WCS parameters in the .qgs file

Finally use this file to view QGIS web client.

I tried to tweak the tool available at https://github.com/oware/Mxd2Qgs. I am successful to certain extent, i.e. creating a .qgs file. But I need some inputs on Getting layer (vector/raster) properties when path for the files is given.

# Import system modules
from xml.dom.minidom import Document
import string
import os
import glob


#Read input parameters from GP dialog

output = "myworld.qgs"

#Create an output qgs file
f = open(output, "w")

# Create the minidom
doc = Document()

# Create the  base element
qgis = doc.createElement("qgis")

qgis.setAttribute("projectname", " ")
qgis.setAttribute("version", "2.4.0-Chugiak")
doc.appendChild(qgis)

# Create the