texts: US Census: bivariate chloropleth of total population (ACS 2013) and total jobs (LEHD 2012)

texts:
The most recent ACS data available is 2013 and Tiger data is also available for 2013. The most recent LEHD data available is 2012 and the geography crosswalk files join it to Tiger 2010 data.

Given that the Tiger files are from different years, is it possible (cartographically sound) to map ACS total population data and LEHD data (WAC, if you're curious) together?


texts: How will Colorblindness impact a GIS career?

texts:
With the wonderful genetics passed from by grandfather via my mother, I'm acquired the blue/ purple-green/brown color blind trait.

How will this affect my capability in properly doing/ interpreting my own or other peoples mapping ?


texts: How to improve accuracy or extrapolate missing or bad GPS points

texts:
Is there an easy to use freeware tool to postprocess my location data, in order to correct
1) location errors (or improve accuracy not based on differential GPS correction) and
2) deal with missing locations, based on ecological information?

e.g. spikes in the location data are indicative of location errors, as are large angles between 3 fixes at small step sizes

Brown hare location data is based on small GPS collar  (U-blox LEA4S), and has been measured every 6 minutes in an open , relatively flat enviroment, including acceleration data every minute for 35 days, 24hr a day


texts: Using VBScript for advanced labeling in ArcMap?

texts:
I am trying to do some advanced labeling with VBScript and need some help getting to the finish line. I have never scripted anything prior to last week so I am definitely a newbie when it comes to VBScripting.

Here is my problem- I have hundreds of soil borings at which several samples have been collected at varying depths. The number of samples collected at each location varies as does the depth. So I wrote a very simple script to get my labels into a table format. I concatenated two fields in excel (Depth and Result) separated by spaces to create the [LABEL] field.

My issue now is adding an expression into my existing script which will return no values if no data is present in the data table. What I'm getting now is just several blank lines when samples were not collected at that depth. See the attached image for reference.

How can I fix my expression in order to skip these blanks?

Here is the script I have thus far...

"<FNT name= 'Arial' size='10'><BOL>" & [Sample_ID] & "</BOL></FNT>" + vbnewline + "<UND>" & "DEPTH (BGS)" & "</UND>" + "  " + "<UND>" & "LEAD (MG/KG)" & "</UND>" + vbnewline + vbnewline + [ALABEL] + vbnewline + [BLABEL] + vbnewline + [CLABEL] + vbnewline + [DLABEL]

This has all been done in the simple label expression maker within ArcGIS so far. I assume I will need to add some kind of if then statement and a loop to get what I need. But like I said earlier... I am no expert at scripting yet.


texts: Object-based Classification in Orfeo ToolBox

texts: I have a multispectral camera with 5 bands (r, g, b, red edge, NIR) and I would like to classify plants object based with Orfeo ToolBox. I am flying with the drone close to the plants and have an image resolution of about 0.01 m per pixel. I guess I have to do a segmentation first, then an object based classification. What steps and models in OTB do I have to follow?


texts: ArcGIS: Geometric Network, Service Areas? OD Cost Matrix?

texts: With ArcGIS Desktop is it possible to run Service areas and OD cost matrix analysis on Geometric networks (ex:gas network)usi ng the network analyst extension?


texts: Error: 'NoneType' object is not iterable (arcpy.ListRasters)

texts:
I am trying to extract a raster through a particular shapefile.

   import arcpy
   from arcpy import env
   from arcpy.sa import *
   arcpy.CheckOutExtension("Spatial")
   env.overwriteOutput = True

   #Set the current workspace
   env.workspace = ("C:\thesis\Forests Only Percent Treecover\2003 BSWM Forest")
   env.nodata = "MINIMUM"
   env.compression = "LZ77"

   #This is for extracting the forest areas of hansen percent treecover
   for raster in arcpy.ListRasters("*0E.tif", "TIFF"):
       print raster #checking the presence of raster
       outputRasterExtractbyMask = ExtractByMask(raster, "NAMRIA.shp")
       outputRasterExtractbyMask_Name = "forests_only"+raster
       outputRasterExtractbyMask.save(outputRasterExtractbyMask_Name)

   print "Finish"

As I print the raster in arcpy.ListRasters to check the presence of raster, this error shows up:

 Traceback (most recent call last):
 File "C:\Users\brentiebark\Dropbox\Python Scripts and mxds\extractbymask.py", line 14, in <module>
 for raster in arcpy.ListRasters("*0E.tif", "TIFF"):
 TypeError: 'NoneType' object is not iterable

I double checked the directory, and my files are really there. Is there something wrong with my code?


texts: Get complete coordinates of a road in OSM Overpass API

texts:
I am trying to get a road entire coordinates using this code

<osm-script output="json">
  <id-query ref="32809713" type="way"/>
  <!-- added by auto repair -->
  <union>
    <item/>
    <recurse type="down"/>
  </union>
  <!-- end of auto repair -->
  <print/>
</osm-script>

The problem with this code is, it only takes part of a road, not entire road. What is wrong?


texts: Not getting different colors on Folium package

texts:
I am trying to plot some points on the map, where the points are marked based on geo-coordinates and each point belongs to a category.
I am trying to set the color of the categories from the color array, and even converted that into a function, yet it is displaying all the categories in one color, as in the image:

Here's the code mentioned below. I am not exactly sure where am I going wrong, since this is from the documentation of Folium.
The variable 'c' here, is actually coming from the column C in dataframe, which has the number of the segment.
mapPlot = folium.Map(location=[latitude, longitude], zoom_start=2)

colors = [
    'red',
    'blue',
    'gray',
    'orange',
    'beige',
    'green',
    'purple',
    'cadetblue',
    'black',
    'pink'
]

def colorCall(val):
  val = int(val)
  return colors[val]

for lon, lat, c in zip(df1['long'], df1['lat'], df1['C']):
    folium.CircleMarker(
        [lat, lon],
        radius=2,

        #fill_color=colors[int(c)],
        fill_color = colorCall(c),
        fill=True,
        fill_opacity=0.7
        ).add_to(mapPlot)

mapPlot


texts:
I have created my own fantasy map, I have the map in a full size png and the map broken down into tiles ( with no geo data )

I installed MS4W and tried to use it to get the georeference data put is seems the program got stuck while trying to do the 4th step here

this is my info I got from the 3rd step:

-Size is 4961, 3231

-Coordinate System is `'

Image Structure Metadata:
-INTERLEAVE=PIXEL
Corner Coordinates:
1st corner Upper Left( 0.0, 0.0 )
2nd corner Lower Left ( 0.0, 3231.0 )
3rd corner Upper Right ( 4961.0,    0.0 )
4th corner Lower Right ( 4961.0, 3231.0 )
Center ( 2480.5, 1615.5 )
Band 1 Block=4961x1 Type=UInt16, ColorInterp=Red
Band 2 Block=4961x1 Type=UInt16, ColorInterp=Green
Band 3 Block=4961x1 Type=UInt16, ColorInterp=Blue
Band 4 Block=4961x1 Type=UInt16, ColorInterp=Alpha

4th step said something about -gcp not available or something

Could someone get my png georeferenced and into a kml or csv file for me so that I can continue my project, or help me with using the right code to continue on

here is my map png


texts:
When I load CityJSON with the CityJSON plugin (QGIS 3.18)
I also see the 2D ground in the 3DMap.
Is that a setting in QGIS?

I can set the 2D symbology to "No Symbols". But this way, I have an empty 2D view.


texts: How to add a TIF raster layer programmatically?

texts:
This is my code whenever i try to load .tif image it shows layer is not valid,QgsRasterLayer always return invalid layer (I am using linux 12.04).

QString myFileName = QFileDialog::getOpenFileName(this,tr("OpenFile"),"/home/ravi",tr("*.tiff"));
std::cout<<"File Name:"<<myFileName.toStdString();
QFileInfo myRasterFileInfo(myFileName);
QgsRasterLayer * mypLayer = new QgsRasterLayer(myRasterFileInfo.filePath(),
myRasterFileInfo.completeBaseName());
if (mypLayer->isValid())
{
   qDebug("Layer is valid");
}
else
{
   qDebug("Layer is NOT valid");
   return;
}


texts: Copying new data into shared Geodatabase in ArcGIS for Server?

texts:
Similar to this question:

Copying new data into shared Geodatabase in ArcGIS for Server using ArcObjects? [closed]

I am wanting to know how to publish a map service without having to copy all of the data to the server because I added one new feature class to our enterprise GDB.

I have gone the route of creating a service definition file without sending up data to the server, and then publishing the sd from Catalog.  However, doing this it seems you have to then delete the existing service (the one I want to replace) and then recreate it from the sd file.  I then have to restart IIS because the indexing in the REST service is off.

Also, with 10.0 we were able to just restart the map service (and IIS) without having to go through all these extra steps and still get the desired result.  Is there a better way to go about this?


texts:
I have some problems with opening a link in ArcGIS API. I have my popup window where LAT, LONG and LINK are shown.

function addPoint(evt) {
        var latitude = evt.mapPoint.getLatitude();
        var longitude = evt.mapPoint.getLongitude();
        var link = "Google";
        map.infoWindow.setTitle("Coordinates");
        map.infoWindow.setContent(
        "lat: " + latitude.toFixed(7) + "
 lon: " + longitude.toFixed(7) + "
 Room: " + link
        );
        map.infoWindow.show(evt.mapPoint, map.getInfoWindowAnchor(evt.screenPoint));

}

Here is the image
However when I click on the link nothing happens (no redirection or opening in new window).
I tried to add window.open() to open the link,

var link = "Google Office";
        on(link, 'click', openLink);
        function openLink(){
            window.open("http://www.google.com");
        }

but then I get the error message

  Uncaught Error: Target must be an event emitter


texts: Extracting attributes for markers in leaflet, onClick event?

texts:
I'm using a JSON file to plot the markers in Leaflet:

 [{
"id": 1,
"date": "1/1/2015",
"Title": "Trinity College",
"Latitude": 41.745167,
"Longitude": -72.69263},
  {
    "id": 2,
    "date": "1/2/2015",
    "Title": "Wesleyan University",
    "Latitude": 41.55709,
    "Longitude": -72.65691
  },{...}]

What I'm doing is the following:

   var markers = new L.markerClusterGroup(); //clustering function

  var markerList = [];

  for (var i = 0; i < jsonDataObject.length; i++) {
    var marker = L.marker(L.latLng(parseFloat(jsonDataObject[i].Latitude), parseFloat(jsonDataObject[i].Longitude)));
    marker.bindPopup(jsonDataObject[i].Title );
    markerList.push(marker);
      }
      markers.addLayers(markerList);
      map.addLayer(markers);

However because I'll add extra functionalities for each individual markers, I want to add 'click' event for each marker and in this function accessing the JSON attributes for each marker. For example:

marker.on('click', onClick_Marker)

function onClick_Marker(e) {
        popup = L.popup()
        .setContent("The number id is: " + e.id)
        .openOn(map);
    }

How to get access individual attributes, from the JSON file, in a pop-up window?


texts: Obtaining extent of each polygon in shapefile using ArcPy?

texts:
In ArcGIS 10 and Python I want to get the extent (xmax, ymax, xmin, ymin) info of each of the polygons in a shapefile.

I can get the extent of the whole shapefile using

file=r"D:\SCRATCH\ARCGIS\100k_trc_tiles_TVM.shp"
desc=arcpy.Describe(file)
print desc.extent.Xmax

  394551.52085039532

But I can't seem to figure out how to get the same info for each row in the dataset.

rows = arcpy.SearchCursor("100k_trc_tiles_TVM")
for row in rows:
 print row

prints the 31 rows in the dataset but

for row in rows:
 desc=arcpy.Describe(row)
 print desc.extent.Xmax

gives an error.

  Runtime error : Object: Describe input
  value is not valid type

I was thinking of adding the extent values to the table using "calculate geometry" but this only gives the centroid. Then I guess we can use something like row.GetValue("xmax").

That being said I know that we can create the X/Y, max/min using the function from http://www.ian-ko.com/free/free_arcgis.htm but it would be best if we can avoid having to add fields, especially if ArcPy can get these values.

Basically I need to get the extents to feed into the clip tool to clip out 30 areas of data (according to the 1:100,000 map sheets) for geoprocessing since the Split tool fails due to the large size of the dataset (see Why Intersect gives ERROR 999999: Error executing function Invalid Topology [Too many lineseg endpoints]?). I want to automate this as it is repeated on a number of datasets.

=== working script ===

# Emulates Arc Info SPLIT tool by using Clip but
# Requires a FC from which each row is used as the input clip feature.
# Each row must be rectangular.
# Used on 12GB FGDB with 100 million records.

#Licence: Creative Commons
#Created by: George Corea; georgec@atgis.com.au, coreagc@gmail.com
import arcpy, string

#inFrame=arcpy.GetParameterAsText(0) # Input dataframe FC
#inFile=arcpy.GetParameterAsText(1) # Input FC for splitting
#outDir=arcpy.GetParameterAsText(2) # Output FGDB

inFrame=r"D:\SCRATCH\ARCGIS\100k_trc_tiles_TVM.shp"
inFile=r"c:\junk\106\data\7_Merge.gdb\FullRez_m2b"
outDir=r"D:\SCRATCH\Projects\206\datasplit\test_slaasp.gdb"
#NameField="Name_1"

#arcpy.env.workspace = r"C:/Workspace"
arcpy.env.overwriteOutput = True

rows = arcpy.SearchCursor(inFrame)
shapeName = arcpy.Describe(inFrame).shapeFieldName
for row in rows:
    feat = row.getValue(shapeName)
    Name = row.Name_1
    print "Executing clip on: "+str(Name)
    extent = feat.extent
    #print extent.XMin,extent.YMin,extent.XMax,extent.YMax
# Create an in_memory polygon
    XMAX = extent.XMax
    XMIN = extent.XMin
    YMAX = extent.YMax
    YMIN = extent.YMin
    pnt1 = arcpy.Point(XMIN, YMIN)
    pnt2 = arcpy.Point(XMIN, YMAX)
    pnt3 = arcpy.Point(XMAX, YMAX)
    pnt4 = arcpy.Point(XMAX, YMIN)
    array = arcpy.Array()
    array.add(pnt1)
    array.add(pnt2)
    array.add(pnt3)
    array.add(pnt4)
    array.add(pnt1)
    polygon = arcpy.Polygon(array)
    ShapeFile = outDir+"\\temp_poly"
    arcpy.CopyFeatures_management(polygon, ShapeFile)

    #print Name
### Set local variables
    in_features = inFile
    clip_features = ShapeFile
    out_feature_class = outDir+"\\"+Name
    xy_tolerance = "0.22"

    # Execute Clip

    try:
        arcpy.Clip_analysis(in_features, clip_features, out_feature_class, xy_tolerance)
        print "Completed: "+str(Name)
    except:
        error = arcpy.GetMessages()
        print "Failed on: "+str(Name)+" due to "+str(error)


texts: Build Interactive Web map for offline/localhost visualisation

texts:
I am part of a small research project where we are leading field experiments in a small watershed, on some remote location far away from our research institute.

I have a series of vector data, shapefiles actually, consisting of various points, lines, polygons, each one representing some specific entity on the field: either rain gage stations, drillings, boreholes, contours, roads, parcels, etc... Also, for those entities, I have some photos taken with a camera to illustrate the actual state of the entity.

I would like to build a series of maps. I want them to be interactive, so that people can choose thematic layer entities of interest and display them, measure distances, etc. Also, it would be good to run them locally on the browser, without the need to install some third party software.

I don't have a clear idea on the way to go. Can anyone help pointing me in the right direction ?


texts: AddIn button disabled in Data View but enabled in Layout View

texts:
VB.NET Addin for ArcMap 10.3 may be causing my Table of Contents to constantly refresh. I have a custom toolbar with a few buttons. I want all of the buttons to be disabled when in Data View and only be enabled in Layout View. I believe that this code works, but is there a better way to do this?

Instead of putting this code on each button, can it be applied to the whole toolbar instead? I put a line of code to write a snippet of text to a log file when the OnUpdate routine is called and the file quickly got huge. That leads me to believe that this code is constantly running.

Public Class btn_DrawBorder
    Inherits ESRI.ArcGIS.Desktop.AddIns.Button
    Private m_Application As IApplication
    Private m_MxDoc As IMxDocument
    Public Sub New()
        m_Application = My.ArcMap.Application
        m_MxDoc = My.ArcMap.Document
    End Sub

    Protected Overrides Sub OnUpdate()
        If TypeOf m_MxDoc.ActiveView Is IPageLayout Then
            Enabled = True
        ElseIf TypeOf m_MxDoc.ActiveView Is IMap Then
            Enabled = False
        End If
    End Sub
End Class


texts: Error when exporting to new folder using ArcPy

texts:
I want to create a backup copy of a feature class and place it in a new folder with today's date written on the folder name. I also need to output file to be a shapefile.

The script below does successfully create the new folder in the style I wanted, but it fails to create the output:

  ExecuteError: ERROR 000210: Cannot create output
  \C:\BackupData\"C:\BackupData\GM_06April2017\LITTER_BINS.shp Failed to
  execute (Select).

import arcpy
import os
import datetime
from arcpy import env

#Set the date to my preferred format ie: 06April2017
Cur_Date = datetime.datetime.now().strftime("%d%B%Y")
print Cur_Date

# Overwrite pre-existing files
arcpy.env.overwriteOutput = True

# Create a new folder called 'GM_' with the current date added to the folder name. Ie. GM_06April2017
newpath = r'C:\BackupData\GM_' + Cur_Date
if not os.path.exists(newpath): os.makedirs(newpath)

#Take a feature class called 'LITTER_BINS' from a File GDB and copy it to the newly created folder, but in shapefile format.
print "Starting backup process..."
env.workspace = r"C:\GIS_folder\Grounds_Maintenance.gdb"
arcpy.Select_analysis("LITTER_BINS", r"C:\Data\GM_" + newpath + "\LITTER_BINS.shp")
print "Bins"

print "Script finished"


texts: Adding location of start and end points to the attribute table

texts: I want to add start and end points to my attribute table for each feature in one layer in terms of actual location and not distance form start meter i.e. start being 0 and end point 10 meters. Is there a way I could do this?


texts: Recursive loop to modify vector layer and refresh mapcanvas

texts:
I'm trying to plot one point "jumping" each second in the canvas (like a gps tracking). I need to take points from a .csv file.

This functionality is made in a QWidget of my application, so the code is the next:

class StartTrack(QWidget):

def __init__(self, MapExplorer):
    QWidget.__init__(self)
    self.MapExplorer = MapExplorer  # MapExplorer is my principal class

def starting(self):                 # I call this method when push a button in the UI
    self.j = 0
    self.addTarget()

def addTarget(self):
    with open('C:\Users\Nacho\SkyLife\GIS4Ship\Archivos\Lex\data\Targets_com.csv') as f:
        self.j += 1
        self.lector = csv.reader(f)
        for i in range(0,self.j):
            row = self.lector.next()
        x = row[1]
        y = row[2]

        # I have this on MapExplorer:

        # self.target_layer = QgsVectorLayer('Point?crs=EPSG:4326', 'points' , 'memory')
        # self.trackprovider = self.target_layer.dataProvider()
        # QgsMapLayerRegistry.instance().addMapLayer(self.target_layer)

        trackprovider = self.MapExplorer.trackprovider
        target = self.MapExplorer.target_layer
        features = []
        feature = QgsFeature()
        feature.setGeometry(QgsGeometry.fromWkt("POINT ("+x+" "+y+")"))
        features.append(feature)
        trackprovider.addFeatures(features)
        target.updateExtents()
        target.triggerRepaint()
        print(self.j, row)

        QTimer.singleShot(1000, self.addTarget)     # Recursive call

My problem is that I can't "delete" the previous point at canvas. I really delete each feature in each iteration, because I'm redefining everything: trackprovider (QgsVectorLayer.dataProvider), target (QgsVectorLayer), features (list of feature), and feature (point).

I want to see only one point each iteration (like moving through canvas) instead of this:

I've already tried to delete feature from list and from layer too, but not results.

Am I missing something related to refreshing canvas, layer or provider?


texts: Getting Error While Parsing Feature To The lineEdit In Python

texts:
am getting error while parsing a selected feature to the lineEdit like, QLineEdit.setText(QString): argument 1 has unexpected type 'long'

layer = self.iface.activeLayer()
    h = layer.selectedFeatures()
    for f in h:
        g =  f['ID']
    self.dlg.lineEdit.setText(g)


texts: Advance Digitizing tools with PostgreSQL views in QGIS

texts:
I have a PostgreSQL database which is accessed by multiple users and each user have different levels of privileges on tables and views. To insert/update/delete features, the users have to work on QGIS.

Let's get to the problem.

I have created a view on MultiLineString data. When i use basic editing tools like add feature, node tool, delete tool and move feature, everything works fine and the edits save successfully. But when i use advance editing tools like split feature, merge features it gives me strange behavior. When i split a linestring and save the edits, one of the splitted segments disappear.  when i merge two segments and save edits, the merged segment disappears.

I tested the advance editing tools i.e. split and merge, on the table from which the above view is created and everything works fine, i am able to save merged and splitted features.

Version: PostgreSQL 9.5 and QGIS 2.18.2


texts: How would I convert NetCDF files to Geotiff format? The data I am using is air temperature data that comes from the Climatic Research Institute and I want to display it on ENVI. I am not familiar with writing codes or scripts.


texts: QGIS - Intersect Spatial GRID and Polygons - Return whole cells on the boundary

texts:
I've been managing within QGIS two shapefiles: the first is a spatial square grid, and the second a city region. Now I would like to intersect the grid keeping only the cells that are overlapping the city, returning the whole cells on the boundaries. So far I used the Intersect function which returns this:

Is there a way (within QGIS or Python, for instance) to keep the whole cells overlapping the boundaries of the blue region?


texts:
I thought I followed the QGIS geo-reference tutorial in precisely.

The following raster (with our overlay) was snapped from GoogleEarthPro (I believe that Google Earth is using EPSG 4326, not WebMercator as is GoogleMaps).

I created 5 place marks, and got the decimal lat/long.

I put them into QGIS 2.10 Georeferencer Picture below.

Tried all variants of settings (Linear, Thin Plate) as well as EPSG 4326, and 3857 output), as well as Linear and Nearest Neighbor Resample.

output is to a GeoTiff.

But the issue still is that you can see that the image is stretched in the X axis.

Now, maybe I should not be viewing in an image viewer (IRFan). And that is my problem.

For when I have (Thin Plate Spline, Linear resample, GeoTiff in EPSG 4326) and view the data in FME Data Inspector, it does look ok!!

All I really want is to read the data out from FME Data Inspector for the Corners points in Lat/Long. So Tell me if I am actually doing things correctly, and the only stupid thing is viewing with IRFAN.

Also given that this is a very small area, what are the best georeference settings for maximum precision?


texts: Functions like iface.actionSelect().trigger() not working in Plugin?

texts:
I want to activate the 'selection Button' in the attribute Toolbar with a click on a button. I always get the error that my class object has no attribute 'iface'.

What am I missing here?

(I am relatively new to building plugins)

I imported everything i need (?):

from PyQt5.QtCore import *
from PyQt5.QtGui import *
from PyQt5.QtWidgets import *
from PyQt5 import uic
from qgis.gui import *
import os
from qgis.core import *
import processing
from qgis.gui import QgisInterface
from qgis.utils import iface

Code looks like this:

self.button.clicked.connect(self.actionSelection)

def actionSelection(self):

    lyr.removeSelection()
    self.iface.setActiveLayer(lyr)
    self.iface.actionSelect().trigger()

lyr.removeSelection works fine, but the two functions below don't!
I get: 'class' object has no attribute 'iface'..

My class looks like this:

WIDGET, BASE = uic.loadUiType(os.path.join(pluginPath, 'Tool.ui'))
class RandomClass(WIDGET, BASE):

I call it like this in my 'main class':

    def run(self):

       self.gui = RandomClass(self.iface.mainWindow())
       self.gui.show()

Could somebody help me?


texts: Create vector .shp file and draw several lines with various attributes (L, Z, S%, etc) from coordinates pairs

texts:
I have a set or coordinates which forms lines from  two points (p1-p2, p2-p3,...., p(n-1)-p(n)), I want to create a SHP file in order to draw several lines from these points, which have several features (Length, Slope, Velocity, etc).

So far, I have been able to create the shape file and create fields, but I can not manage to fill these fields or create a geometry lines to associate with.
Can anybody point me to an example  to successfully draw lines on a shapefile using a foor loop, because coordinates points  are in the order of magnitude of 10^2.

#IMPORTAR LIBRERIAS
import ogr,os
import osgeo.osr as osr

#DRIVERNAME
driver = ogr.GetDriverByName('ESRI Shapefile')
ds = driver.CreateDataSource("PROYECTO_" + shp_name + os.path.sep + "00_GIS" + os.path.sep + "00_IN" + os.path.sep + shp_name + '.shp')

#REFERENCIA ESPACIAL
srs = osr.SpatialReference()
srs.ImportFromEPSG(32717) #WGS84 17S

#CREAR SHAPE de PUNTOS (ogr.wbkPoint) O LINEA (ogr.wkbLineString)
if shp_type == "linea":
    layer = ds.CreateLayer(shp_name, srs, ogr.wkbLineString)
if shp_type == "punto":
    layer = ds.CreateLayer(shp_name, srs, ogr.wkbPoint)

#AGREGAR CAMPOS
#NOMBRE TRAMO (Numero o Caracter)
fieldDefn = ogr.FieldDefn('Ramal', ogr.OFTString)
layer.CreateField(fieldDefn)

#NOMBRE TIPO DE ALCANTARILLADO (Sanitario, Pluvial, Combinado)
fieldDefn = ogr.FieldDefn('Tipo', ogr.OFTString)
layer.CreateField(fieldDefn)

# NOMBRE TIPO DE ALCANTARILLADO (Sanitario, Pluvial, Combinado)
fieldDefn = ogr.FieldDefn('Material', ogr.OFTString)
layer.CreateField(fieldDefn)

# LONGITUD
fieldDefn = ogr.FieldDefn('Longitud', ogr.OFTString)
layer.CreateField(fieldDefn)


texts:
If we have coded domains in the geodatabase model, what are any disadvantages to changing the codes to the Description. For example:

Existing Coded Domain:

Code, Description
1, Single Family Dwelling
2, Multi Family Dwelling
3, No Tenants

OR

Code, Description
5200, Concrete
5201, Steel
5175, Copper

New Coded Domain with Code and Description matching:

Code, Description
Single Family Dwelling, Single Family Dwelling
Multi Family Dwelling, Multi Family Dwelling
No Tenants, No Tenants

OR

Code, Description
Concrete, Concrete
Steel, Steel
Copper, Copper


texts: Merging OSM files with OSMOSIS and import to Postgres with OSM2PGSQL

texts:
I have a problem merging and then importing two OSM files.
I downloaded to distinct .osm-files from cloudmade. They are not sharing any geographical attribute (as the first is the country Monaco, the 2nd is the county Bremen in Germany).

I merged the with OSMOSIS and it gives me an output file.

Then I import it via OSM2PGSQL.

Looking at the dataset using QGIS it either shows me the data of Bremen or Monaco (depending the order in the OSMOSIS merge statement).

Does anyone has a solution for this problem?

Thanks a lot
Carsten


texts: How to reduce transparency for Openlayers using Map Composer?

texts: Increasing transparency of Openlayers will increase the contrast of others layers. For sure, it will increase the esthetical aspect of maps.


texts: Multiple ring buffer in ArcMap gives Update parameters execution error

texts:
My issue is with multiple ring buffer in ArcGIS Desktop 10.8.1. I'm dealing with GIS1 workbook Tutorial 9-3. Every time trying to enter distances I get this error message telling me:

Update parameters execution error`

What should I do?


texts: How to calculate polygon centroids in R (for non-contiguous shapes)

texts:
I've spent a little while figuring out the answer to this question.
It's not immediately obvious from a Google search, so thought it may useful to post the answer on here. There is also an additional question about non-contiguous polygons.

Instant easy answer: use the command:

centroids <- getSpPPolygonsLabptSlots(polys)

(This was found in the class description of the SpatialPolygonsDataFrame R data class for the overarching spatial package in R, sp)

This seems to do exactly the same thing as

cents <- SpatialPointsDataFrame(coords=cents, data=sids@data, proj4string=CRS("+proj=longlat +ellps=clrk66"))

in the following code, which should be replicable on any R installation (try it!)

#Rcentroids
install.packages("GISTools")
library(GISTools)
sids <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1],
                      proj4string=CRS("+proj=longlat +ellps=clrk66"))
class(sids)
plot(sids)
writeSpatialShape(sids, "sids")
cents <- coordinates(sids)
cents <- SpatialPointsDataFrame(coords=cents, data=sids@data,
                  proj4string=CRS("+proj=longlat +ellps=clrk66"))
points(cents, col = "Blue")
writeSpatialShape(cents, "cents")

centroids <- getSpPPolygonsLabptSlots(sids)
points(centroids, pch = 3, col = "Red")

Where cents (blue) and centroids (red) are identical centroids (this should plot should appear after you've run the code):

So far so good. But when you calculate polygon centroids in QGIS (menu: Vector | Geometry | Polygon Centroids ), there are slightly different results for non-contiguous polygons:

So this question is 3-things:

A quick and easy answer
A warning for people using R to calculate centroids for non-contiguous polygons
A question about how it should be done in R to properly account for multi-part (non-contiguous) polygons


texts: QGIS feature count is wrong after changing categories programmatically

texts:
I want to delete categories in my vector-layer style if the feature-count for the category is zero.

The following code works so far if you run it in the python console but the feature count for the layer is wrong after the categories with zero objects have been deleted.

If i turn on the edit mode, create a feature in the layer and close the attribute-input-dialog the feature count is recalculated.

How can I trigger the refresh of the feature count programatically?

# encoding: utf-8

from PyQt4.QtGui import *
from PyQt4.QtCore import *
from qgis.core import QgsVectorLayer, QgsFeature, QgsFeatureRequest, QgsExpression, QgsCategorizedSymbolRendererV2
from qgis.gui import QgsMessageBar

layer =iface.activeLayer()
anzahl=0
deleted=0
if type(layer) == QgsVectorLayer and layer.hasGeometryType():
    r = layer.rendererV2()
    if r.type() == 'RuleRenderer':
        anzahl = len(r.rootRule().children())
        # code snippet from g.carillo: https://gis.stackexchange.com/questions/209674/qgis-test-rules-with-python
        for rule in r.rootRule().children(): # Iterate rules
            if rule.filter(): # Filter out empty rules
                request = QgsFeatureRequest( rule.filter() )
                count = len( [f for f in layer.getFeatures( request )] )

                if count <1:
                    r.rootRule().removeChild(rule)
                    deleted +=1
    elif r.type() == 'categorizedSymbol':
        anzahl = len(r.categories())
        counter=0
        del_categories=[]
        new_categories=[]
        for position, category in enumerate(r.categories()): # Iterate rules
            myVariable = category.value()
            expression_string = u" \""+unicode(r.legendClassificationAttribute())+"\" = '"+myVariable+"' "
            #print expression_string
            filter_string = QgsExpression( expression_string)
            request = QgsFeatureRequest( filter_string )
            count = len( [f for f in layer.getFeatures( request )] )

            if count <1:
                del_categories.append(category)
                deleted +=1
            else:
                new_categories.append(category)
                print layer.name() + ' : ' + category.label() + ' : ' + myVariable + ' : ' + str(count)

            counter+=1
        print del_categories

        r.deleteAllCategories()

        for i in new_categories:
            try:
                r.addCategory(i)
            except:
                print "oh no"

    # apply the renderer to the layer
    if r is not None:
        layer.setRendererV2(r)
    layer.dataProvider().forceReload()
    layer.triggerRepaint()

    iface.layerTreeView().refreshLayerSymbology(layer.id())

    print("Layer style filtered! "+str(deleted)+" of "+str(anzahl)+" rules filtered")


texts: Optimization of multi-band raster clip by polygon

texts:
I have tiled (50x50), 365 band, rasters in my postgis database. Each band represents a day in a year.
I want to return the raster pixels which intersect with an input geometry, including all bands. I'd like to figure out the most optimized way of doing this.
Current Approach:

filter the raster table with ST_Intersects(rast,geom)
Clip raster to geom with ST_Clip(rast,geom)
St_Dumpvalues(clip) the result of clip to examine the band-pixel values of the clip

Function:
CREATE OR REPLACE FUNCTION get_time_series (start_year INT, end_year INT, geom geometry)
RETURNS TABLE (
    year INT,
    day_of_year INT,
    result float8[]
)
AS $$
BEGIN
RETURN QUERY
 with myresult as(
        with myclip as (
            SELECT rast_table.rid,st_clip(rast_table.rast, geom) as clip,rast_table.year
                FROM rast_table_50x50 rast_table
            WHERE ST_Intersects(rast_table.rast, geom)
                AND rast_table.year BETWEEN start_year AND end_year
        )
        SELECT (st_dumpvalues(clip)).*,myclip.year from myclip
    )
    SELECT myresult.year,myresult.nband as day_of_year,valarray as result FROM myresult;

END; $$
LANGUAGE 'plpgsql';

Is there another approach which might yield better peformance?
Would a plpython function be of any use for optimizing the clip/mask operation?
Should I be clipping at all?


texts: Overlaying WMS from Mapserver in OpenLayers with epsg:900913 from Google/Bing/OpenStreetMap

texts:
ADDED: Another curious effect is that no matter to what scale I zoom on the baselayer google/bing map, the WMS overlay is always the same-sized shot of MN.  That is, as the baselayer map is scaled up or down, the WMS overlay map never scales -- it always shows the same bounds.

I am trying to convert an Openlayers page from running overlay maps with OpenLayers.Layer.MapServer to one running WMS with OpenLayers.Layer.WMS.  The information I am providing through these interfaces is overlaying on a Google/Bing/Openlayers baselayer.  The Mapserver implementaion works fine, but the WMS implementation seems to generate an overlay from my data that is ~ 2x too large for the overlay with the baselayer maps.

** Please excuse the use of links here and for pastebin... I wanted to embed images and make links directly to pastebin, but being a new user, I was given restrictions of no images and only two links.

Here is a picture of the working overlay:

Here is a picture with the problematic WMS overlay:

The mapfiles for mapserver are mostly identical, as I used the working Mapserver mapfile to create the mapfile for WMS, but with some necessary mods to the WMS mapfile to make it WMS enabled.

Here is the working Mapserver map file:
http://pastebin.com/VQBW5YSx

And here is the WMS mapfile:
http://pastebin.com/swMagqma

Any help is appreciated.

Thanks,
Matt


texts: I’m trying to make a point file with elevations using the tutorials I’ve found on youtube. The conversion seems to work like it should but when I import the points into a CAD program the contours are correct but they are separated an equal distance vertically. Is there a setting I’m missing? Can post a picture if need be.


texts: Is it possible to change the legend so the symbol is on the right side?

texts: Is it possible to change the legend arrangement so the symbol will be in right side and the label description will be in the left side? The QGIS default is label description in the right side and the symbol is in the left side


texts: Register an image using mutual information

texts:
I want a procedure for registering two images based on their mutual information. I took a pair of georeferenced radar images of the same location, edited the metadata of one so its position was off by about a km, then tried to use mutual information in the ENVI image registration workflow to correct the position of the edited image so it matched the unedited image. The result was a horrible warp that did not correct the problem at all.
The images have near identical conditions (there's a very slight change in brightness values, but nothing that should be a problem).

How can I use mutual information in ENVI to register a pair of images?


texts: Why are spatial lagged exogeneous variables rarely used?

texts:
In the publication "Regional Productivity and Income Convergence in the  Unified Germany 1992-2000" (link) it is claimed without detailed explanation: "Usually,  in  spatial  econometric  models  only  spatial  lags  in  the dependent variable are taken into account, while spatial lagged exogenous variables ordinarily are not explicitly modelled." And:  "As  in  time  series  analysis  one  can  argue  that  spatial  effects  stemming  from  exogenous variables will be captured by a spatial lagged endogenous variable."

From these statements I get that both approaches are in principal equal. Is this true? If it is, why exactly is the first approach (spatially lagged endogenous/dependent variable) preferable in the spatial domain?


texts:
I have a CHM model (raster) with values -10 up 100.

I need to reclassify values as:

value < 0 => 0
0 < value < 50 => the same value
value > 50 => 0

I tried so many tools, but I cannot find the correct way. Do you have any ideas?


texts: Adding data to a feature layer on ArcGIS Online using arcrest and arcresthelper.featureservicetools?

texts:
I'm attempting to automate some web map updates on ArcGIS Online through the use of the arcrest module.  The following two methods are intended to act as the interface between local data sources and AGOL.  The first of the two methods, "TruncateWebLayer" works as intended and removes all of the existing records from the specified feature layer.  The second method, "AddFeaturesToWebLayer" always results in an error.

#This always works correctly
def TruncateWebLayer(agolID = '', lyrName = ''):
    securityinfo = {}
    securityinfo['security_type'] = 'Portal'
    securityinfo['username'] = 'username'
    securityinfo['password'] = 'password'
    securityinfo['org_url'] = 'https://myportal.maps.arcgis.com'
    securityinfo['proxy_url'] = None
    securityinfo['proxy_port'] = None
    securityinfo['referer_url'] = None
    securityinfo['token_url'] = None
    securityinfo['certificatefile'] = None
    securityinfo['keyfile'] = None
    securityinfo['client_id'] = None
    securityinfo['secret_id'] = None

    fst = featureservicetools.featureservicetools(securityinfo)
    fs = fst.GetFeatureService(itemId=agolID,returnURLOnly=False)
    fl = fst.GetLayerFromFeatureService(fs=fs,layerName=lyrName,returnURLOnly=False)

    fst.DeleteFeaturesFromFeatureLayer(fl.url, '1=1')

#This always fails on the AddFeatures line
def AddFeaturesToWebLayer(addFeatures, agolID = '', lyrName = ''):
    securityinfo = {}
    securityinfo['security_type'] = 'Portal'
    securityinfo['username'] = 'username'
    securityinfo['password'] = 'password'
    securityinfo['org_url'] = 'https://myportal.maps.arcgis.com'
    securityinfo['proxy_url'] = None
    securityinfo['proxy_port'] = None
    securityinfo['referer_url'] = None
    securityinfo['token_url'] = None
    securityinfo['certificatefile'] = None
    securityinfo['keyfile'] = None
    securityinfo['client_id'] = None
    securityinfo['secret_id'] = None

    fst = featureservicetools.featureservicetools(securityinfo)
    fs = fst.GetFeatureService(itemId=agolID,returnURLOnly=False)
    fl = fst.GetLayerFromFeatureService(fs=fs,layerName=lyrName,returnURLOnly=False)

    #This always fails..
    fl.addFeatures(fc=addFeatures, attachmentTable=None)

    #I've tried this too and it fails with the same error message because it references addFeatures also.
    #fst.AddFeaturesToFeatureLayer(r'url to feature layer', addFeatures, chunksize=0, lowerCaseFieldNames=False)

The error message I always receive is as follows:

Traceback (most recent call last):
  File "<string>", line 181, in <module>
  File "<string>", line 74, in AddFeaturesToWebLayer
  File "C:\xxxxx\xxxxx\xxxxx\xxxxx\xxxx\xxxx\site-packages\arcrest\agol\services.py", line 2330, in addFeatures
    for chunk in chunks:
  File "C:\xxxxx\xxxxx\xxxxx\xxxxx\xxxx\xxxx\site-packages\arcrest\agol\services.py", line 1972, in _chunks
    l.sort()
TypeError: '<' not supported between instances of 'dict' and 'dict'

This looks like the error has something to do with sorting of the chunk dictionaries in the the services.py module.  I cannot make sense of this.  Is this maybe not the correct method for doing what I'm intending it to do?  I have noticed that running these methods from an ArcGIS Pro session and replacing the addFeatures line with an arcpy.append_management works and everything gets uploaded to the web layer, however, it takes an very long time when dealing with feature layers with thousands of records.


texts: I am 13 and I have an assessment on G.I.S soon and I don't have a clue what it is.


texts: How can I close the identify result dialog programmatically?

texts: I am trying to close the identify result dialog which appeared when I select to identify a feature from the active layer. Does someone know how can I do it?


texts: PostGIS Multipolygon when converted to shapefile doesn't match the source geometry

texts:
I am trying to convert a MULTIPOLYGON record from PostGIS (v2.2.2) to the shapefile via QGIS (v2.18.9), but the source geometry and the resulting shapefile do not match as displayed on the images below (the red layer is the source and the green layer is the resulting shapefile). I get the same result if I perform the conversion via ogr2ogr from GDAL v2.0.0. This conversion takes around 40 minutes to complete.

If I dump the source geometry to POLYGON and perform the conversion, the resulting shapefile is correct and the conversion is much faster (7 seconds + 17 seconds for dump). But, instead of 5 features, now I end up with 305188.

The source MULTIPOLYGON records are created by ST_Union on the adjacent cells if that matters.

Am I missing something? Is there a way to perform the correct conversion? Can anyone explain why is this happening?

Source geometry

Exported shapefile with filled holes marked by arrows

I have more details:

pgsql2shp produces the same resulting shapefile like QGIS or ogr2ogr.

@dbaston - I did cut of most of the data and left the part around the area presented on the above images (around ('POINT(-89.45 29.99)',4326)). I converted this small sample and got the same result. If you or anybody wants to try the conversion on the same sample, you can get the dump here:

Sample data backup from PostgreSQL v9.5

After 31 hours, the ST_IsValid returned true on all geometry from the complete data set. I got the same result for the smaller sample available on the link above. However, Check validity in QGIS for the same smaller data set returned:

ring 1 of polygon 712 not in exterior ring
ring 2 of polygon 712 not in exterior ring
ring 3 of polygon 712 not in exterior...

How come there is a difference in the validity check in PostGIS and QGIS?

To check further, I removed all the polygon parts except the one marked blue on the image above. I checked its WKT format and it looks fine. And now the conversion to shapefile with only this feature was correct.

Does this make any sense? I did check more datasets created by ST_Union on the adjacent cells and found out that the same issue exists there too when converting to shapefile. The same data converted to GeoJSON produces the correct result.

MULTIPOLYGON was created from POLYGON geometry with the following expression:

st_multi(st_union(st_buffer(geom,0)))


texts: ReportLab integration into QGIS returns error message "no module named reportlab"

texts:
I want to use ReportLab in my own plugin which I have written for QGIS. I have installed ReportLab for the python versions 2.4, 2.5 and 2.6. The package is respectively installed  into the site-packages folder (e.g. C:/Python24/Lib/site-packages/reportLab).

If I want to use the tool in my plugin, the plugin returns the error "ImportError: no module named reportLab". I have set the path to the modul into the environment variables, but this doesn´t help.

I try to import the module by "from reportLab.pdfgen import *".

Do you have any other suggestions or ideas?


texts: Find centroids of rows with Multipoint Geometry

texts:
I wish to find the centroid of a set of given points. The points will be rows of my SpatiaLite database and have Multipoint geometry. I checked the SpatiaLite manual but the ST_Centroid function I find says it takes as input a Surface s (while the corresponding function in PostGis takes as input Multipoint as well).
I thought of creating a surface from the Multipoint rows but it seems I have to do it in this order: Multipoint->Linestring (with MakeLine)->Area (BuildArea)->Centroid (with ST_Centroid), something that seems a heavy workaround to even try. Of course I tried it but it does not work.
Can anyone suggest a fast workaround?


texts:
I am attempting (and failing) to add an image to a map layout using ArcObjects via Python comtypes. I do need to run this outside of ArcMap (i.e. not in the ArcMap Python window).

The basic idea is: open mxd, place jpg in layout, save/close mxd.

Python solution preferred, but I would also accept a similar workflow in VB.NET.

Using ArcGIS 10.2.2, Python 2.7.5

Script:

from comtypes.client import GetModule, CreateObject
from Snippets import GetLibPath, InitStandalone

def CType(obj, interface):
   """Casts obj to interface and returns comtypes POINTER or None"""
   try:
       newobj = obj.QueryInterface(interface)
       return newobj
   except:
       return None

def main():
    InitStandalone()
    modCarto = GetModule(r"C:\Program Files (x86)\ArcGIS\Desktop10.2\com\esriCarto.olb")
    modGeom = GetModule(r"C:\Program Files (x86)\ArcGIS\Desktop10.2\com\esriGeometry.olb")
    modArcMapUI = GetModule(r"C:\Program Files (x86)\ArcGIS\Desktop10.2\com\esriArcMapUI.olb")
    modFramework = GetModule(r"C:\Program Files (x86)\ArcGIS\Desktop10.2\com\esriFramework.olb")

    app = CreateObject(modFramework.AppROT, interface=modFramework.IAppROT)
    mxDoc = CType(app.Item(0).Document,modArcMapUI.IMxDocument)
    layout = mxDoc.PageLayout
    graphicContainer = CType(layout, modCarto.IGraphicsContainer)

    pictureElement = CreateObject(modCarto.JpgPictureElement, interface=modCarto.IPictureElement)

    envelope = CreateObject(modGeom.Envelope, interface=modGeom.IEnvelope)
    envelope.PutCoords(0, 0, 5, 5)

    pictureElement.Geometry = envelope
    pictureElement.ImportPictureFromFile(r'C:\junk\happy.jpg')

    iElement = CType(pictureElement,modCarto.IElement)

    graphicContainer.AddElement(pictureElement,0) # ERROR HERE
    #graphicContainer.AddElement(iElement,0)

    app.Item(0).SaveDocument()

if __name__ == '__main__':
    main()

Traceback (whether attempting to add pictureElement or iElement):

Traceback (most recent call last):
  File "C:\Scripts\arcobjects_add_jpeg.py", line 54, in <module>
    main()
  File "C:\Scripts\arcobjects_add_jpeg.py", line 47, in main
    graphicContainer.AddElement(pictureElement,0)
COMError: (-2147024809, 'The parameter is incorrect.', (None, None, None, 0, None))

Full disclosure: I am absolutely new with ArcObjects, so whether or not I am handling/creating/juggling objects between interfaces correctly, I am not sure, but I think I'm close.


texts:
I am attempting to allow for the pop-up to be toggled on and off by pressing a checkbox. Here are some snippets of the code:

var clickHandler

function popupOn(evt){
var query = new esri.tasks.Query();
query.geometry = pointToExtent(map,evt.mapPoint,10);
var deferred = featureLayer.selectFeatures(query,esri.layers.FeatureLayer.SELECTION_NEW);
map.infoWindow.setFeatures([deferred]);
map.infoWindow.show(evt.mapPoint);
}

function clickConnect(connect){
if(connect){
//perform the identify task on click
clickHandler = dojo.connect(map, "onClick", popupOn);
}
else{
//disconnect the click handler
dojo.disconnect(clickHandler);
clickHandler = null;
}
}

<input type='checkbox' data-dojo-type="dijit.form.CheckBox" checked="checked" data-dojo-props="onClick: dojo.connect('true')"/>Tick here to turn selecting on.


texts: Splitting *.osm.pbf file into multiple files?

texts:
I need to split a PBF file into multiple PBF files.

Is this possible with any tool using Python?

I looked at osmosis for example, but from what I understand from the documentation it is only possible to split the file based on bounding boxes/polygons. I would need to split it into relatively equal parts without knowing any underlying geography. The multiple files should all be valid.

The destination of the data should be other PBF files. The initial PBF can be any geographic area (mostly continents or world). The separate PBF files would then be processed in parallel on different machines. As such I am not really bothered about geographic boundaries and it would be nice not to worry about them. I just need to generate a bunch of valid PBF files. Valid should also means that there should be no information lost from the original file.


texts: PDAL pipeline could not connect to server: Connection refused

texts:
I am trying to import .las files to the postgres database. My pipeline is:

{
  "pipeline":[
{
  "type":"readers.las",
  "filename":"data/pt000002.las",
  "spatialreference":"EPSG:4326"
},
{
  "type":"filters.chipper",
  "capacity":500
},
{
  "type":"writers.pgpointcloud",
  "connection":"host='localhost' dbname='db_example' user='postgres' password='super'",
  "table":"las2",
  "compression":"dimensional",
  "srid":"4326"
}
  ]
}

This pipeline is worked for me before when I have PDAL installed in Conda environment. Now, I have PDAL installed in Docker environment, and when I call pipeline with:

sudo docker run -v /home/user/Documents/projectX:/data pdal/pdal:1.7 pdal pipeline /data/laz2pg.txt

I get the error:

 pdal/pdal:1.7 pdal pipeline /data/laz2pg.txt
 PDAL: could not connect to server: Connection refused
 Is the server running on host "localhost" (127.0.0.1) and accepting
 TCP/IP connections on port 5432?
 could not connect to server: Address not available
 Is the server running on host "localhost" (::1) and accepting
 TCP/IP connections on port 5432?

All answers are true. Anyway, I used the same database when I used PDAL with conda. I am using Ubuntu 18.04 Bionic.


texts: Extract source points in r.cost.full for each gridcell

texts:
R.cost.full from the GRASS algorithm calculates an accumulated cost surface (i.e. a map of travel costs) for a unit cost layer/friction map. The output is a raster map where each grid cell has a value that represents the costs for getting to the closest source (e.g. the next town). Is there a possibility to extract which town is the closest source of each grid cell?

My intention is to get a second raster that contains the closest source (town) for all grid cells.

The manual states that by setting the parameter:
outdir=name I can create a raster to contain movement directions. However this is not what I exactly need.


texts: Copy and paste a non-spatial record in the same table (make a duplicate a record)

texts:
How can I duplicate (copy/paste) a non-spatial record into the same table in ArcMap? I can copy and paste to/from spatial tables, but there doesn't seem to be a way to do this in the attribute table of a non-spatial table (there's no paste button).

I need a way to copy and paste, as a common editing task (multiple times per minute). Maybe with one or two clicks, and no noticeable processing time. I would have thought I'd be able to copy and paste, just as I do with many other computer programs.

I am using ArcGIS Desktop 10.3.1.


texts: Geoserver / SLD multiline lables not working?

texts:
for several hours Im trying to create a multiline label for a line.
  I use Geoserver 2.5, the line is in postGIS, the style is very simple :

<?xml version="1.0" encoding="UTF-8"?> <sld:StyledLayerDescriptor

xmlns="http://www.opengis.net/sld" xmlns:sld="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc" xmlns:gml="http://www.opengis.net/gml" version="1.0.0">   <sld:NamedLayer>
        <sld:Name>Borders_1fe760e8</sld:Name>
        <sld:UserStyle>
          <sld:Name>Borders_1fe760e8</sld:Name>
          <sld:Title>Borders</sld:Title>
          <sld:IsDefault>1</sld:IsDefault>
          <sld:FeatureTypeStyle>
            <sld:Name>name</sld:Name>
            <sld:Rule>
              <sld:Name>rule1</sld:Name>
              <sld:Title>Title 1</sld:Title>
              <sld:Abstract>-</sld:Abstract>
              <sld:LineSymbolizer>
                <sld:Stroke>
                  <sld:CssParameter name="stroke">#FF0000</sld:CssParameter>
                  <sld:CssParameter name="stroke-width">2</sld:CssParameter>
                </sld:Stroke>
              </sld:LineSymbolizer>
              <sld:TextSymbolizer>

                  <ogc:PropertyName>Name_1</ogc:PropertyName><![CDATA[
]]><ogc:PropertyName>Name_2</ogc:PropertyName>

                <sld:Font>
                  <sld:CssParameter name="font-family">Tahoma</sld:CssParameter>
                  <sld:CssParameter name="font-size">16</sld:CssParameter>
                  <sld:CssParameter name="font-style">normal</sld:CssParameter>
                  <sld:CssParameter name="font-weight">normal</sld:CssParameter>
                </sld:Font>
                <sld:LabelPlacement>
                  <sld:PointPlacement>
                    <sld:AnchorPoint>
                      <sld:AnchorPointX>0.0</sld:AnchorPointX>
                      <sld:AnchorPointY>0.5</sld:AnchorPointY>
                    </sld:AnchorPoint>
                    </sld:PointPlacement>
                </sld:LabelPlacement>
                <sld:Fill>
                  <sld:CssParameter name="fill">#FF0000</sld:CssParameter>
                </sld:Fill>
                <sld:VendorOption name="followLine">true</sld:VendorOption>
              </sld:TextSymbolizer>
            </sld:Rule>
          </sld:FeatureTypeStyle>
        </sld:UserStyle>   </sld:NamedLayer> </sld:StyledLayerDescriptor>

yet I am unable to see the lables 1 and 2 on two lines - they always end up on the same line.
  I tried :

SLD - CDATA with line break (as seen above)
SLD - CDATA with
 tag
SLD - use only one field, POSTGIS field value with \n (value1\nvalue2)
SLD - use only one field, POSTGIS field value with
 tag (value1
value2)

and many other permutations of above, i copy-pasted the code, typed it, used geoserver style editor, typed it manually into xml file directly on the server,  but I wasnt able to produce two line label. I did try all the tricks I was able to find via google.

What am I doing wrong ? Or is it some known Geoserver bug ?


texts: Get Database Table Name From Layer Name in ArcObjects

texts:
I am using ArcObjects v10.0 with a SQL Server 2008 database.  Given the name of a Layer, is there a way for me to extract the corresponding name of the database table that represents that Layer?

The Layer name and the Database Table name will often be identical, but they will not in cases where the Layer name contains invalid characters, such as spaces.


texts: Spatialite queries not loading entirely on canvas

texts:
I connected a SQLite database to my QGIS project. I wrote a query in the DB manager whose result contains 261 elements.

Then, I tried to load the result on the canvas (from the load button at the bottom of the query window) but nothing happened.

So, I created a view (which still contains 261 elements) which appeared in the table list of the database. Then, I added the view to the canvas but only 100 elements are loaded. I don't understand what happened and how to fix this.

Does anyone have an idea?


texts: Unable to create topology with pgrouting

texts:
I am just trying to create a sample topology using pgrouting.  The dataset I am using is simply two crossing line segments (not connected in the middle because in the future I will have over/underpasses and I want to see how this is handled).

However I keep getting this error:

NOTICE:  PROCESSING:
NOTICE:  pgr_createTopology('test_edges', 0.1, 'id', 'geom', 'source', 'target', rows_where := 'true', clean := f)
NOTICE:  Performing checks, please wait .....
NOTICE:  ----> PGR ERROR in pgr_createTopology: Column source not found
HINT:    ----> Check your column name
CONTEXT:  SQL statement "SELECT _pgr_onError(err, reportErrs, fnName,  'Column '|| col ||' not found', ' Check your column name','Column '|| col || ' found')"
PL/pgSQL function _pgr_getcolumnname(text,text,text,integer,text) line 18 at PERFORM
SQL statement "select *                 from _pgr_getColumnName(sname, tname,source,2,fnName)"
PL/pgSQL function pgr_createtopology(text,double precision,text,text,text,text,text,boolean) line 61 at SQL statement
NOTICE:  Unexpected error raise_exception
Total query runtime: 64 msec
1 row retrieved.

The table name is "test_edges" and it has two columns, "id" and "geom".  I am using this query:

select pgr_createTopology('test_edges', 0.1, 'id', 'geom');

This is the table I am using:

Does anybody know why this would not work?


texts: Create a center node based on both polygon and line shapefile

texts: Sorry for cross posting. I want to develop a node layer from two other layers. One of the layers is a polygon contains a series of zones. Another layer is a line layer of roads. The nodes I want to create are basically, for each zone, centers of the roads instead of center of zone. For example, if a zone contains a park in the north with no roads there (zone 2), the node should be on the south side based on the roadway density in zone 2. Am I able to do that?


texts: How to rubberband zoom w/o the shift key: arcgis api for javascript

texts:
I want to enable the rubberband zoom effect without users having to hold the shift key down.  This behavior is found is Esri navigation bar dijit.  I don't want to use the dijit.

How can I get a zoom in/out bounding box to draw without forcing users to use shift key?


texts: Generating Service Areas in Network Analyst with Added Roads?

texts:
Here’s what I am running into with Network Analyst- the goal is to create service areas for fire stations based on travel times.

•   I have a baseline county wide roads shapefile that I have used as input to successfully create a network dataset and create subsequent service areas based on travel times.

•   Next, I add new roads (i.e. digitize by hand) to the baseline roads dataset for a part of the county where a new housing development is going in.  I am then able to create/build a network dataset without any hitches.  However, when I try to create a new service area using the updated network dataset based on the updated roads shapefile, the new roads that I have added are not included/taken into account in the service area polygons that are generated.  I've entered valid values for mph, length (in miles), and travel time (in minutes) for the new records.

•   In short, my goal is to be able to add new roads to an existing roads shapefile and have those roads included when network analyst solves for my service area query.

•   I have run the Integrate tool on the shapefile with the new roads and I still receive the same results - the new roads are not incorporated in the service area polygons that are produced.

I am fairly new to using network analyst.


texts: Canopy density & NDVI Regression using ArcGIS for Desktop?

texts:
I have a Landsat ETM+ image, I also have 200 sample points with coordinates,canopy density, different band values and NDVI values.

If I build OLS model I only got sample points regressed.

How can I join the image and the points so I can do regression and produce a map in ArcGIS for Desktop?


texts:
I created in OL3 selector elements new.ol.interacion.Select () on a WFS layer, and when I try to take these features returns an empty array.

My intention is to create a json with the selection check the map features.

This is my code:

  var Urbana_source = new ol.source.Vector({
   format: new ol.format.GeoJSON(),
   url: function(extent, resolution, projection) {
       return 'http://localhost:8080/geoserver/wfs?service=WFS&' +
    'version=1.1.0&request=GetFeature&typename=pucol:Urbana&' +
    'outputFormat=application/json&srsname=EPSG:3857&' +
    'bbox=' + extent.join(',') + ',EPSG:3857';
      },
     strategy: ol.loadingstrategy.tile(ol.tilegrid.createXYZ({
     maxZoom: 19
     }))

     });

 var urbana_wfs = new ol.layer.Vector({
 source: Urbana_source,
 style: new ol.style.Style({

  stroke: new ol.style.Stroke({
   color: 'rgba(0, 0, 0, 0)',
   width: 2
   })
 })
});
 urbana_wfs.set('name', 'urbana');

select = new ol.interaction.Select({
        condition: ol.events.condition.click
    });
map.addInteraction(select);
//var feature_buff = select.getFeatures();
select.getFeatures().on('change:length', function(e)
{
    var feature_buff = select.getFeatures();

    var geojson  = new ol.format.GeoJSON();
    var  json    = geojson.writeFeatures(feature_buff);

    });

and image


texts:
When needing to represent True / False or Yes / No which of the following examples is better from a geodatabase design point of view?

Our database architect hates the option on the right but when looking at the data using a non-esri tool using letters makes it easy to understand without needing to do a database join.

Store an integer      or   Store a string and use the first letter

0  - 'No' or                        'N' - 'No'

1  - 'Yes' or                       'Y' - 'Yes'

0  - 'False' or                       'F' - 'False'

1  - 'True' or                   'T' - 'True'


texts:
If I try, like this way:

        $.ajax({
       jsonp: false,
       jsonpCallback: 'getJson',
       type: 'GET',
       url: "http://79.141.150.167/geoserver/xxx/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=xxx:business_service&maxFeatures=50&outputFormat=text/javascript&format_options=callback:getJson",
       async: false,
       dataType: 'jsonp',
       success: function(data){
          console.log(data);
       }
    });

My browser tells me following:

Resource interpreted as Script but transferred with MIME type text/xml: "http://79.141.150.167/geoserver/xxx/ows?service=WFS&version=1.0.0&requ…tputFormat=text/javascript&format_options=callback:getJson&_=1393412236410".

Uncaught SyntaxError: Unexpected token < (ows:1)

...what I'm I doing wrong?


texts: Error with Annotate Selected Features in ArcGIS Pro 2.1.2

texts:
I'm trying complete a script that will automatically create an annotation layer from a feature layer in ArcGIS Pro but am hitting a wall with the "Annotate Selected Features" code block (all information is on the ESRI site here).

Here's my code:

import arcpy
import os
import os.path
import sys
import datetime

arcpy.env.workspace = r"C:\Users\myname\Desktop\USACE\Working_Map\Working_map.gdb"
aprx = arcpy.mp.ArcGISProject(r"C:\Users\myname\Desktop\USACE\Working_Map\Working_Map.aprx")
m = aprx.listMaps("Map")[0]
lyr = m.listLayers("WGRES_Layer")[0]

arcpy.AnnotateSelectedFeatures_cartography(m, 'WGRES_Layer', "WGRES_LayerAnno 'Class 1'", 'GENERATE_UNPLACED')

Now the weirdest thing is the error message I get when I try to run this block; I get the message:

Traceback (most recent call last):
   File "<string>", line 15, in <module>
AttributeError: module 'arcpy' has no attribute 'AnnotateSelectedFeatures_cartography'

Does anyone have any idea what's going on?

I get the same thing if I tried to run the tool as "AnnotateSelectedFeatures" instead (i.e. remove "_cartography" from the end).


texts: Importing xyz text ascii GPS time-series for velocity visualization using QGIS?


texts: Strange "RuntimeError: ERROR 160195: An invalid SQL statement was used". when trying to debug SearchCursor

texts:
I'm trying to get a SearchCursor to function in my script. According to everything I could find, this should have worked. The code reads like this:
with arcpy.SearchCursor(sickTrees,'OID@') as cursor:
    for numb in distList:
        print('hi')

Edit: This is the full code now
import arcpy
arcpy.env.overwriteOutput = True
sickTrees=r"C:\Users\Owen Hansen\Desktop\Python GIS\Ch5\Chapter5_LabData\Chapter5_assignmentData.gdb\sicktrees"#arcpy.GetParameterAsText(0)
fineTrees=r"C:\Users\Owen Hansen\Desktop\Python GIS\Ch5\Chapter5_LabData\Chapter5_assignmentData.gdb\healthytrees"#arcpy.GetParameterAsText(1)
distList=[5,10,15,20,25]
test='OID@'
with arcpy.da.SearchCursor(sickTrees,test) as cursor:#It is on fine, help, this bloody SQL error won't go away
for numb in distList:
    with arcpy.da.SearchCursor(sickTrees,test) as cursor:#fixed it by moving the cursor creation
        count=0#should reset, isn't. Only cycling once
        print(str(numb))
        for tree in cursor:
            print('hi'+" "+str(count)+' '+str(tree[0]))
            count+=1
    del cursor
test=['OID@']
with arcpy.da.SearchCursor(sickTrees,test) as cursor:#tried for loop, was dumb
    for numb in distList:
        percentAvgCount=0.0
        for tree in cursor:#tried buffer, was dumb
            #select by oid tree[0] could select by attribute
            arcpy.management.SelectLayerByAttribute(sickTrees,'NEW_SELECTION','OBJECTID = '+str(tree[0]))
            arcpy.management.SelectLayerByLocation(sickTrees, 'WITHIN_A_DISTANCE', sickTrees, numb, 'NEW_SELECTION', 'NOT_INVERT')#honestly, I'm just guessing how this actually works
            sickTreeCount=float(arcpy.management.GetCount(sickTrees).getOutput(0))
            #select by oid tree[0]
            arcpy.management.SelectLayerByAttribute(sickTrees,'NEW_SELECTION','OBJECTID = '+str(tree[0]))
            arcpy.management.SelectLayerByLocation(healthyTrees, 'WITHIN_A_DISTANCE', sickTrees, numb, 'NEW_SELECTION', 'NOT_INVERT')
            helTreeCount=float(arcpy.management.GetCount(healthyTrees).getOutput(0))
            treePercent=sickTreeCount/(helTreeCount+sickTreeCount)
            percentAvgCount+=treePercent
    totalAvgPercent=percentAvgCount/float(numb)
    print(str(numb)+' %:'+' '+str(totalAvgPercent))
del cursor


texts: Plotting X, Y coordinates generated from MySQL database in QGIS

texts:
Currently running QGIS 1.8.0 on Mac OS X 10.8.4.

I can plot points on a map using the "Add delimited text layer" plugin to get x-y coordinates from a text file. However, I'd like to be able to extract the x-y coordinates from a MySQL database directly rather than having to create an intermediate text file. Can this be done? The database consists of numerous tables containing info about postcodes, eastings, northings, and other information. I can select "Add vector layer" and connect to a MySQL database but I'm not given the option to provide the SQL query necessary to extract the coordinate data or to define which variables should be defined as x and y coordinates.


texts: Creating Selection MapTool in PyQgis Custom Application

texts:
I'm developing a custom GIS application using PyQgis API.

How can I let the user select features on the map using Maptools? I would like to replicate the select feature there in Qgis to make the selection on click of a button. My understanding is to use layer.selectedFeatures() function , but I am not sure of how to add the Selection Maptool to the mapcanvas and toolbar.


texts:
Is there a dataset for the spatial extents of historic levels of the Great Lakes since the Last Glacial Maximum? More specifically, we are looking for the maximum level of the Great Lakes after 5,000 years before present.

I'm collaborating on some research on distribution of an aquatic species in the Great Lakes area. I am trying to see if the current distribution can be explained by the historic water levels of the Great Lakes.

I thought I would check to see if anyone is aware of such a dataset before I look into creating something myself.


texts: Change of GPS firmware Geoexplorer Trimble

texts:
Is it possible to install the firmware of a GPS of another brand, the Geoxplorer XH 2005 of Trimble?

Informative note: the Geoxplorer XH 2005 has Windows Mobile Operating System 5.1


texts: Shortest route between large amount of points

texts:
I have to check the street lights of various communities. There are approx. 540 lamps to check, I was wondering if QGIS (or any other opensource software like PostGIS or R) can help me finding the best (aka shortest) way. I have the lights-to-check in a shapefile.
I guess it is too complex to calculate with my PC but thought it might be worth a shot to ask here for any creative solutions!


texts: Making multiple legend boxes all same size using ArcGIS for Desktop?

texts: How do I make all the rectangle or (legend box) all the same ? Look at the circles I marked in the pictures. It should be very simple thing to do. Anyone out there have a good workaround ?


texts: Programmatically Load Hyperlink File To Event Layer in ArcMap 10.2

texts: I am using an arcpy script tool to create an XY Event Layer from a text file and then save the layer to a .lyr. This works, but I also want to load a hyperlink file and apply it to the layer automatically. I tried using UpdateLayer, but this just replaces the entire layer, it doesn't keep the same source data. So, does anyone know if there is a way to load a .hlk file programmatically? The solution does not have to be in Python.


texts: How to call a GP service from ArcObjects in .NET using ArcGIS Engine 10.2?

texts:
I published and successfully accessed a GP service in ArcMap 10.2, but I couldn't make it work in my .NET program using ArcObjects 10.2:

        IGeoProcessor gp = new GeoProcessorClass();
        gp.AddToolbox(@"http://localhost:6080/arcgis/rest/services;mygpfolder/testsel");

        //input
        IGPUtilities2 gputils = new GPUtilitiesClass();
        IFeatureClass fc = (IFeatureClass)gputils.OpenDatasetFromLocation(@"E:\input.shp");
        IRecordSetInit recordset = new RecordSetClass();
        recordset.SetSourceTable((ITable)fc, null);
        IGPRecordSet gprecordset = new GPFeatureRecordSetLayerClass();
        gprecordset.RecordSet = (IRecordSet)recordset;

        IGPFeatureRecordSetLayer recordlayer = new GPFeatureRecordSetLayerClass();
        recordlayer = gprecordset as IGPFeatureRecordSetLayer;

        IVariantArray parameters = new VarArrayClass();
        parameters.Add(recordlayer);

        gp.Execute("Select", parameters, null);

I got "Error HRESULT E_FAIL has been returned from a call to a COM component." at "gp.Execute".  The GP service is published using tool "Select" under toolbox "Analysis" in ArcMap:

Any suggestions?


texts:
I need to update fields in attribute table of a vector using QGIS. I use field calculator for each field, but i have too many fields (and too many vectors also) and i need to program the task. The code that i'm using for each first field is:

  segALD00*(((xmin($geometry) - xmax($geometry))^2+(ymin($geometry) -
  ymax($geometry))^2)^0.5)*102.47/longkm

to the last field

  segALD23*(((xmin($geometry) - xmax($geometry))^2+(ymin($geometry) -
  ymax($geometry))^2)^0.5)*102.47/longkm

Then another vector:

  segETH00*(((xmin($geometry) - xmax($geometry))^2+(ymin($geometry) -
  ymax($geometry))^2)^0.5)*102.47/longkm

to the last field

  segETH23*(((xmin($geometry) - xmax($geometry))^2+(ymin($geometry) -
  ymax($geometry))^2)^0.5)*102.47/longkm

How can i code this?


texts: Duplicating attributes table within QGIS composer?

texts:
I'm working with QGIS 2.14.6-Essen on a Windows 7 OS.

I have added an attributes table in QGIS composer with specific properties (attributes displayed, fonts, frame and background colors, etc.). This table filters data from a layer ("code" = '086077').

Now, I want to add a second attributes table with exactly the same specifications except the filter ("code" = '086233').

Is there a way to duplicate the first table so that I would just have to modifiy the filter value in the duplicated table? So far, I have no other choice than starting from the beginning for my second table...


texts: Determining if road is privately owned or not in USA

texts:
I'm working on a project where it matters greatly whether or not a piece of land is on a private road.

Google Maps in many cases takes private roads (for example in master-planned subdivisions) and treats them in the same way it does for public roads.

Is there any quick way to check whether a road is privately owned or public?

In a lot of cases, even if a road is public it can be privately owned.


texts: I am working on layering GIS over time axile to visualize vehicle's movement in time

texts: I am trying to combine time and latitude/longitude dimensions by layering calendar over map and essentially creating time-space map showing how vehicle moves in time and over the roads. Anybody have any experience, what kind of application I can use?


texts:
I'm trying to do a spatial join on a PostGIS table that has a bunch of hexes and each hex has an ID. I have a separate table withs points in lats and longs and I want to know which hex each point belong in. For instance (points are made up)

 id     lat    long
 1     -83.1    24.5
 2      55.9   -90.2
 3      16.7    23.4

to

id      lat     long   hexid
 1     -83.1    24.5     5
 2      55.9   -90.2    87
 3      16.7    23.4    23

The problem I have is that I have alot of points (over 100 million) and doing the spatial join is taking DAYS. I have broken the query up into special geographic regions and running them concurrently, but is there anything I can do to my query that will improve my performance? My query is below:

CREATE TABLE newtable AS select id, lat::double precision, long::double precision

ALTER TABLE newtable ADD Primary key (id);

SELECT AddGeometryColumn ('public','newtable','geom',4326,'POINT',2);

UPDATE newtable SET geom = ST_SetSRID(ST_MakePoint(long, lat) ,4326);

ALTER TABLE newtable ADD COLUMN hexid numeric;

UPDATE newtable set hexid = hextable.hexid from hextable WHERE  ST_Within(newtable.geom, hextable.geom);


texts:
I am using EasyButton plugin of leaflet (link). More specifically I applied EasyBar which unites several buttons. At certain zoom levels, I would like the buttons to collapse, in other words displaying only one button.

I have seen examples but couldn't succeed to apply in mine.
Here is my code:
So, how can I "hide" buttons (3,4,5,6) except one (2) in certain zoom levels?

 var btn2 = L.easyButton ({position: 'topright',
  states: [{
    stateName: 'NoSources2',
    icon:'masquer le nombre de resources',
      title: 'NoSources' ,
      onClick: function(control) {
                map.removeLayer(NoSources);
                control.state('NoSources1')},
    },
    {
     stateName: 'NoSources1',
    icon: 'voir le nombre de resources',
    title: 'NoSources',
    onClick: function(control)  {
    map.addLayer(NoSources);
    control.state('NoSources2')},
  }]
});
btn2.addTo(map);

var btn3 = L.easyButton ({position: 'topright',
  states: [{
    stateName: 'altitude2',
    icon:'masquer l`altitude des chapelles',
      title: 'altitude' ,
      onClick: function(control) {
                map.removeLayer(altitude);
                control.state('altitude1')},
    },
    {
        stateName: 'altitude1',
    icon: 'voir l`altitude des chapelles',
    title: 'altitude',
    onClick: function(control)  {
    map.addLayer(altitude);
    control.state('altitude2')},
  }]
});
btn3.addTo(map);

var btn4 = L.easyButton ({position: 'topright',
  states: [{
    stateName: 'orientation1',
    icon:'masquer l`orientation des chapelles',
      title: 'orientation1' ,
      onClick: function(control) {
                map.removeLayer(orientation);
                control.state('orientation2')},
    },
    {
        stateName: 'orientation2',
    icon: 'voir l`orientation des chapelles',
    title: 'orientation1',
    onClick: function(control)  {
    map.addLayer(orientation);
    control.state('orientation1')},
  }]
});
btn4.addTo(map);

var btn5 = L.easyButton ({position: 'topright',
  states: [{
    stateName: 'cyclique1',
    icon:'masquer l`info des temps cyclique',
      title: 'cyclique' ,
      onClick: function(control) {
                map.removeLayer(cyclicTime);
                control.state('cyclique2')},
    },
    {
        stateName: 'cyclique2',
    icon: 'voir l`info des temps cyclique',
    title: 'cyclique',
    onClick: function(control)  {
    map.addLayer(cyclicTime);
    control.state('cyclique1')},
  }]
});
btn5.addTo(map);

var btn6 = L.easyButton ({position: 'topright',
  states: [{
    stateName: 'linear1',
    icon:'masquer l`info des temps linéaire',
      title: 'linear' ,
      onClick: function(control) {
                map.removeLayer(linearTime);
                control.state('linear2')},
    },
    {
        stateName: 'linear2',
    icon: 'voir l`info des temps linéaire',
    title: 'linear',
    onClick: function(control)  {
    map.addLayer(linearTime);
    control.state('linear1')},
  }]
});
btn6.addTo(map);

var buttons = [btn4, btn3, btn2, btn6, btn5];
// build a toolbar with them
L.easyBar(buttons, {position: 'topright'}).addTo(map);


texts: Break down geoprocessing operation (Buffer) into smaller subsets

texts:
I have a feature class with 790,000 points (and it is constantly increasing). I need to select some of them (approx. 350,000), and then buffer that selection.

Currently, the process runs like this (pseudocode):

arcpy.MakeFeatureLayer_management(fcIn, tempLayer, '"CUST" LIKE \'%G%\'')
arcpy.Buffer_analysis(tempLayer, fcOut, "150 Feet", "FULL", "ROUND", "ALL", "")

The Buffer step is prone to failure due to memory errors, which isn't terribly surprising. To get around that, I would like to further subset the selection, creating a bunch of buffers and then merging them all at the end. I know that this will take at least as long as (and probably longer than) a single Buffer operation, but avoiding a complete failure is a higher priority.

The question: how do I structure a loop that iterates selecting 10,000 (or so) points at a time and then buffers them? I am struggling to turn the iteration concept into code that will work with feature layers.

(And of course, if this is a terrible, unworkable idea, please let me know. I'm not very committed to it yet.)


texts: How can I use GK + ETRS89 projected data in QGIS or Proj?

texts:
I got a dataset of Hamburg, Germany in so called "Lagestatus 320" (LS320). The metadata describes it as "Gauss-Krueger mit ETRS89 [Lagestatus 320]". Coordinates are like 3565575.3,5936627.0, they are metric.

The epsg database does not have such a coordinate reference system. The closest option is EPSG:31467 aka "DHDN / 3-degree Gauss-Kruger zone 3". Of course, this uses the wrong ellipsoid.

If I assign that to my data, things end up ~600 meters north of their proper location.

How can I properly use this data in QGIS or any other proj-based software?


texts: How to add a rule to topology in ArcMap?

texts:
Looks like I have the same problem as Anthony, but it seems like I have a bit more information to it.

--> Running ArcGIS 10.4.1

import os
import arcpy
arcpy.env.overwriteOutput = 'TRUE'
gdb = r'F:\path\to\gdb\database.gdb'
fd = 'check_geometries'
fc = 'geometries'
topo = 'geometries_topology'
arcpy.env.workspace = gdb
arcpy.CreateTopology_management(fd, topo)
topo_path = os.path.join(gdb, fd, topo)
fc_path = os.path.join(gdb, fd, fc)
arcpy.ValidateTopology_management(topo_path)
arcpy.AddRuleToTopology_management(topo_path, "Must Not Intersect (Line)", fc_path)

The error code I am getting here is continuously 000732.

File "C:\APPS\ArcGIS\Desktop10.4\ArcPy\arcpy\management.py", line 16214, in AddRuleToTopology
raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Input Feature class: Dataset   does not exist or is not supported
Failed to execute (AddRuleToTopology).

Strange enough, the same dataset is used to create the topology in the first place. I am using the name of the feature class to create the topology and the complete path of the topology to add the rule because it was failing before, when only giving the name, and I read that it requires the full path.

Would also be great if ESRI would get that consistant...

In addition to the earlier referenced posting I added the ValidateTopology function, which runs through without any problem. Still, adding a rule to the topology is impossible.

I can see the error is raised in line 16214 in the management.py file. Since I went step by step through the process I was able to confirm that all paths went correct through the gp_fixargs() function, and optional argument not given have been also correctly taken of the arguments list.

I assume the process is failing in the subsequent run of gp.AddRuleToTopology_management().

--> AddRuleToTopology

16208   from arcpy.geoprocessing._base import gp, gp_fixargs
16209   from arcpy.arcobjects.arcobjectconversion import convertArcObjectToPythonObject
16210   try:
16211       retval = convertArcObjectToPythonObject(gp.AddRuleToTopology_management(*gp_fixargs((in_topology, rule_type, in_featureclass, subtype, in_featureclass2, subtype2), True)))
16212       return retval
16213   except Exception, e:
16214       raise e

Can someone tell me on how to get a rule added to my topology?


texts:
How to connect a text in the website to a polygon on leaflet map?
(I am not refering to labeling - a simple html text)
Or a more advanced option: how to highlight this polygon when mouse is over the text?


texts: How to move features from one layer to another?

texts: I have a polygon which belongs to a layer and I need to move it to another one. Is there a way to do it without drawing it again?


texts:
I wrote a script that, ideally, loads raw .las (LiDAR data) files into a geodatabase, then creates a terrain from the multipoint files that are created.  When I run the script though I get:

  ExecuteError: ERROR 000210: Cannot create output C:\Users\peter\Desktop\test.gdb\bare_earth\1las
  The table name is invalid.
  Failed to execute (LASToMultipoint)

So far the only thing I can figure out is that Windows 7 keeps changing my folder permissions to read-only, but I can still manually create, edit, and delete in the directories.

import os
import arcpy, math, glob
from arcpy import env

arcpy.CheckOutExtension("3D")
arcpy.env.overwriteOutput = True
#make workspace in memory
arcpy.env.workspace = "in_memory"

#define parameters using GetParameter in ArcMap
geodatabase = arcpy.GetParameter(0)         #0)input Geodatabase
outfeaturedataset = arcpy.GetParameter(1)   #1)Desired Name for the feature class
projection = arcpy.GetParameter(2)          #2)coordinate system
terrain_name = arcpy.GetParameter(3)        #3)Name of the Terrain
class_code = arcpy.GetParameter(4)          #4) the return class code to use
lasdirectory = arcpy.GetParameter(5)        #5) directory to get las data in

#designate directory where .las files are located
os.chdir (lasdirectory)

created_FC_tool = arcpy.CreateFeatureDataset_management(geodatabase, outfeaturedataset, projection)
created_FC = created_FC_tool.getOutput(0)
#Create a Terrain

arcpy.CreateTerrain_3d (created_FC, terrain_name, 1.24605815948, 50000, '', 'WINDOWSIZE', 'ZMIN', 'NONE', 1)

output_terrain = os.path.join (created_FC, terrain_name)
arcpy.AddTerrainPyramidLevel_3d (output_terrain, 'WINDOWSIZE', '2.5 1000; 5 2000; 10 4000; 20 8000; 40 16000')
#process the .las file and add to the feature dataset and terrain

for LASfiles in glob.glob("*.las"):
    print LASfiles
    pointinfo = "point_file_information"
    #Get Pointfile information
    arcpy.PointFileInformation_3d(LASfiles, pointinfo, "LAS", "", projection)

    #get average point spacing variable
    pointspacing = arcpy.SearchCursor (pointinfo)
    field = "Pt_Spacing"
    for row in pointspacing:
        avg_space = (row.getValue(field))
    print avg_space

    #convert las  to multipoint

    #create lable for multipoint file
    listlable = list(LASfiles)
    del listlable[-4]
    multipointlable = ''.join(listlable)
    print multipointlable

    print created_FC
    outputPlable = os.path.join (created_FC, multipointlable)
    arcpy.LASToMultipoint_3d(LASfiles, outputPlable, 1.5, class_code, "", "INTENSITY", projection, "las", 1)

    print outputPlable

    outputPlable_parameters = outputPlable + ' Shape masspoints 1 0 40 true false R06305328las_embed <None> false'
    print outputPlable_parameters
    arcpy.AddFeatureClassToTerrain_3d (output_terrain, outputPlable_parameters)

#Build the terrain
arcpy.BuildTerrain_3d (output_terrain)

If it is a Windows 7 problem I don't know what I'm going to do as my work computer is locked down by our help desk.


texts: Oracle spatial: Point cluster vs MultiPoint?

texts:
I'm reading Oracle's documentation and i can't see the difference between Point cluster and MultiPoint element type.

The following example code is from the documentation

INSERT INTO t1 (i, d, g)
VALUES (
  16,
  'Point cluster',
  sdo_geometry (2005, null, null, sdo_elem_info_array (1,1,3),
    sdo_ordinate_array (50,5, 55,7, 60,5))
);
INSERT INTO t1 (i, d, g)
VALUES (
  17,
  'Multipoint',
  sdo_geometry (2005, null, null, sdo_elem_info_array (1,1,1, 3,1,1, 5,1,1),
    sdo_ordinate_array (65,5, 70,7, 75,5))
);

It seems that:

Point cluster is one element (of many points -in that example 3 points: 50,5, 55,7, 60,5) as it needs one triplet in sdo_elem_info_array

MultiPoint is many elements ( and each element is a Point -in that example 3 points: 65,5, 70,7, 75,5) as it needs three triplets in sdo_elem_info_array

And further down documentation refers to them as the same thing(?):

2.5.8 Several Geometry Types
Example 2-13 creates a table and inserts various geometries,
including multipoints (point clusters), multipolygons, and collections...

What is this mess? What's the difference? The sdo_geometry requires different sdo_elem_info_array parameters so there must be a difference between those types?


texts: How to distinguish SRTM from DTED using C#?

texts: Does anyone know how to easily distinguish between SRTM vs DTED programmatically without opening file metadata header? both are .dt1 or .dt2 extensions.


texts: Applying Zoom Extent from several QgsVector layers with Labels in PyQGIS?

texts:
I have the following case:

Project consists of three layers:

OpenLayers plugin with Google streets
QGSVector layer with Point features and Label enabled
QGSVector layer with Polyline features and Label enabled

I need to (1) zoom extent canvas for both #2 and #3 and (2) make sure Labels are visible for all the features placed on these two layers.

I'm using the following code:

        idList = []
        idList.append(vlPoint.id())
        idList.append(vlLine.id())
        render = QgsMapRenderer()
        render.setLayerSet(idList)
        canvas = qgis.utils.iface.mapCanvas()
        canvas.setExtent(render.fullExtent())
        canvas.refresh()

The result of this code execution is below:

As you see the very upper point feature is displayed partially and its label is missed at all.  If I zoom out canvas manually just a little- everything is displayed OK.

Need to know how to correctly zoom to extent in my case - to make sure all the features and labels are displayed.

If I have only one vector layer - the following code solves such issue:

 vlPoint.selectAll()
 canvas.zoomToSelected()
 vlPoint.removeSelection()

But for two layers this approach doesn't work.


texts: GPS Module Without Velocity or Altitude Limitations

texts:
I am looking for a GPS module that I can use in a high-altitude project.  The problem I have discovered is that GPS manufacturers are required to place a set of restrictions on GPS modules.  The restrictions set a maximum speed and altitude at which the device should stop functioning.  Basically, this limits someone's ability to build a cruise missile.

From what I understand, it is not illegal to have a GPS receiver without these limitations, it is illegal to export it from the U.S. without considering it a weapon.

If I am not mistaken, GPS is also used to track some satellites in orbit, so as long as the orbiting satellite is not in geo-synchronous orbit or further.  Obviously, these GPS receivers do not have the limit, as both the height and velocity would far exceed the mandate.

Does anyone know of a GPS unit that does not have these limitations, or, is there a unit that allows you to flash the firmware and manually remove the limitations?

Link to mention of limitations: http://en.wikipedia.org/wiki/CoCom#Legacy


texts: Adding a store to GeoServer: Restricting access to file system

texts:
When adding a store to Geoserver I can restrict access to the underlying file system by modifying web.xml (ref: https://docs.geoserver.org/stable/en/user/security/webadmin/filebrowse.html). However, is it possible to restrict access for a given user in GeoServer to a specific directory on the underlying server?

I've looked into Geofence, but as far as I can tell all users added to GeoServer will still be able to browse to all of the GeoServer data when adding a store. Preferably I would like to specify access to a particular directory on the server for a given GeoServer user when they add a new store. Is this possible?

Thanks in advance for any suggestions.

I'm using GeoServer 2.16.0 on Ubuntu 18.04.


texts: How to query all records (>1000) and export to csv?

texts:
I posted in arcgis forum but never got any answer.
Can anyone tell me how to solve this?

Assume you have a vehicle tracking records from Oct.1 to Dec. 31 with wind speed information.
All data is saved in geodatabase (sde-sqlserver) and up in ArcGIS Server 10.1 as feature layer.
A time slider show the locations of a car with wind speed.

When a user change a time-range (eg. Oct.2 - Oct.4), first query task (for count) calculate the number of features within the range.
There are usually more than 1000 results even two days (eg. 1750) (I don't want to change this limit though).

I used another query task (executeforIds) to keep all records in hand but reduce the amount with modulus (1/10) which is still enough number for making a nice chart for an overall trend of wind speed.
However, I also want to provide an option to download a whole dataset in csv (in this case 1750 rows)

Here, I used findtask to retrieve attribute dataset within the time-range.

try {
    //console.log(app.objIDs);
    var attribs;
    csvdata = "";

    var find = new esri.tasks.FindTask("http://xyz.com/arcgis/rest/services/Realtime/Car_Wind_WM/MapServer");
    var params = new esri.tasks.FindParameters();
    params.layerIds = [0];
    params.searchFields = ["OBJECTID"];

    for (var i = 0; i < app.objIDs.length; i++) {
        params.searchText = app.objIDs[i];
        find.execute(params, function (results) {
            attribs = results[0].feature.attributes;
            csvdata += attribs.Date_Central + ',' + attribs.Longitude + "," + attribs.Latitude + "," + attribs.windspeed + "\n";
            console.log(csvdata); //1
        }, function (error) {
            alert("Error");
        });
    }
} catch (error) {
    alert("Change the time range first");
}
console.log(csvdata); //2
setTimeout(function () {
    formatData(app.csvdata);
}, (3 * 1000));

csvdata was undefined at console.log2, although all csvdata (in console.log1) show data in console.
The issue seems timing of find.execute complete, so I added setTimeout.

This seems work but when i increase the time-range, it won't obviously.

Is there anyway to keep all records (1000 - 200,000) in a specific time range and export to csv?


texts: Allowing identification tool to click on symbol instead of feature in QGIS

texts:
Working on a rather crowded map, I would like to ease the identification of colocated point features, each one linked to a photo. I usually have 2 or 3 photos on the same point and would like to keep this granularity as at some point I will add to each photo its source geometry and heading.
What I would like to achieve is (below a specific map scale threshold) to explode those points and let the user click on each photo by clicking on its symbol instead of the feature itself.
Otherwise the user has to click on the point, then click on the correct feature in the identification results to make sure it's the correct one, then right-click + open formular. 4 clicks instead of one is quite cumbersome and can take quite some time if this action is frequent.
This is one use case but in general I think it would be useful to identify a feature by simply clicking its symbol. There was already a close use-case in a feature request long ago, but I am wondering if there might be some easier solution first.
Any ideas?


texts:
I need to use https://github.com/kartena/Proj4Leaflet for my case.
I have EPSG:32638
But I don't know what parameter to pass to constructor. From the web site
it has example:
var crs = new L.Proj.CRS('EPSG:2400',
  '+lon_0=15.808277777799999 +lat_0=0.0 +k=1.0 +x_0=1500000.0 ' +
  '+y_0=0.0 +proj=tmerc +ellps=bessel +units=m ' +
  '+towgs84=414.1,41.3,603.1,-0.855,2.141,-7.023,0 +no_defs',
  {
    resolutions: [8192, 4096, 2048] // 3 example zoom level resolutions
  }
);

I am confused about the second and third parameters, what should they be?
Also once I pass this CRS to my map, can it mess up my map? (It has already many features which I implemented with Leaflet)


texts: Listing Elements or Datasets of Layer gives ERROR 000840 value is not Feature Layer?

texts:
I used Ar5cPy to do batch processing of Network Analyst - Make Service Area layer of thousand of stores.

Everything went well except now, I am stuck with layer (.lyr) files and I would like to conduct further analysis on the datasets that they contain.

Basically I would have

"SA_Store_001.lyr"  (and so on..)

And when I open this inside ArcMap, it has the usual elements of the Service Area Analysis.

including

"Facilities"
"Route Barriers"
"Polygons"
etc.

Suppose, I would like to SelectAttributes by polygons, can you give me a skeleton code for this? The help files on ArcGIS on has SelectLayerByAttribute -- what I want is the Polygons.

Unfortunately, this didn't work:

Even the first line doesn't work (because I don't know how to deal with the .lyr file type)

inputLayer = 'D:\\ARCGIS\\FoodSystemsAnalysis\\Coverage\\Auckland\\SMC_SM_Ta9.lyr'

arcpy.MakeFeatureLayer_management(inputLayer, "SMC_Ta9")

  arcgisscripting.ExecuteError: Failed to execute. Parameters are not
  valid. ERROR 000840: The value is not a Feature Layer. Failed to
  execute (MakeFeatureLayer).

I think the MakeFeatureLayer_management takes features as input.


texts: PermissionError when trying to load the layers in an ESRI File Geodatabase

texts:
I am trying to list the layers of a file geodatabase using fiona. I would like to then pick a layer and use geopandas to extract the layer as a geodataframe:
import fiona

layers = fiona.listlayers(r'main_datasets/BRT/TOP10NL.gdb')

I get an error as follows:

fiona.errors.DriverError: Permission denied

I did the following:

Removed any .lock files which exist within the gdb
Tried opening other gdbs to see if I got the same error

Removing lock files in this particular gdb did not help. Strangely, the method in fiona successfully lists layers in other gdbs which don't have any notable differences (I checked in ArcCatalog and these are also FileGDB's as far as I can tell).
Why is this particular gdb giving me trouble? What could be causing the permission error?


texts: Show raw bytes that make up an attribute in a shapefile

texts:
I am processing shapefiles, and I'm trying to track down an encoding problem. I want to check the raw bytes of an attribute to see if this shapefile (and hence the step that created it) is the source of the problem.

If I have a shapefile, how can I print out the raw bytes for attribute X for the object which has attribute Y = Z? e.g. I suspect that the row ID=12 has a 'bad' encoding for the NAME field. How can I print out the raw bytes (e.g. a hexidecimal format or od -c format, or ...) that make up the NAME field for the object in the shapefile with ID=12? If it's only possible to print out all attributes for that object that's OK (There aren't many and I can visually search for what I want).

Preferably with open source tools / things installable on ubuntu, and preferably on the command line.


texts:
I am currently working with a shapefile made up of multi-part polygons symbolizing zoning areas. Geocortex is a web-based service used to display your data. When I click on a larger polygon in ArcMap the polygon highlights neatly along the lines as a selected feature should. When I do the same through the Geocortex web-based data viewer, the polygon selected highlights from random points - creating triangular connection points - and not the line represented. I'll show an example in the attached image.

If anyone has experienced this issue, then how can I fix it?


texts: How to do a spatial join in QGIS with two layers with different CRS?

texts:
I am trying to do a Spatial Join (Link Attributes based on location) with two shapefiles. What I want to do is determine which streets are in a certain region, or in other words are contained in that region.

I have a drawn the shapefile of regions which will be the join layer, but I have realised too late that this shapefile was not in the project CRS which is Belge lambert 72. The shapefile I have drawn as the join layer is WGS 84.

I have CRS transformation on so QGIS is actually projecting it in the right place it should be.

However, when I have tried to to change the CRS from WGS 84 to lambert 72, either it does not work, or the location changes?

I thought it would not change anything but I was wrong, as you can see in the table above all value's that are joined to the street layer are NULL.

How do I solve this?