texts: ArcPy Feature Class to Feature Class blank results - ArcGIS Pro

texts:
I have this Feature Class to Feature Class script that creates a feature class of patient visits for each store for every month for every year it has been open. However, when I run the script it will create feature classes by month and year with the correct labels, for example (Jan2015_NWPatientVisits2000Gresham), but when I look at the attribute table it contains no data.

I'm sure I am missing something small, but I just can't seem to figure out what it is.

Here is my script:

import arcpy

# Set environment settings
arcpy.env.workspace = r"C:\arcGIS_Shared\Python\CenterHeatMaps.gdb"

#Declare variables
fc = 'Open_Store_Centers'
fields = ['USER_market_id','USER_Store_ID','USER_Store_Center_Name']
fieldname = 'USER_market_id'

#Define WHERE clause statement
whereclause = """{} = 2000""".format(arcpy.AddFieldDelimiters(fc, fieldname))
sqlclause = (None, 'Order By USER_market_id, USER_Store_ID')

# loop through months
years = [2015, 2016, 2017]
months = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']

for year in years:
    for month in months:
        with arcpy.da.SearchCursor(in_table = fc, field_names = fields, where_clause=whereclause, sql_clause=(None, 'ORDER BY USER_market_id, USER_Store_ID')) as cursor:
        #Loop through each row established in cursor
            for row in (cursor):
            # Set local variables for FeatureClasstoFeatureClass
                inFeatures = "PatientVisitsGeocoded"
                outLocation = r"C:\arcGIS_Shared\Python\CenterHeatMaps.gdb"
                outFeatureClass = "{2}{3}_NWPatientVisits{0}{1}".format(row[0], row[2], month, year)
                delimitedfield = arcpy.AddFieldDelimiters(arcpy.env.workspace,"USER_CenterID")
                expression = delimitedfield + "= {0}".format(row[1]) + " AND USER_DOSYear = {0}".format(year) + " AND USER_DOSMonth = '{0}'".format(month)

#         Execute FeatureClassToFeatureClass
                arcpy.FeatureClassToFeatureClass_conversion(inFeatures, outLocation, outFeatureClass, expression)

        #Print Results
                print(row[2])
                count = arcpy.GetMessageCount()
                print (arcpy.GetMessage(count-1))


texts: Using if/elif/else statements in ArcGIS Field Calculator?

texts:
I have no experience with Field Calculator in ArcGIS 10.

I have a very simple calculation I would like to perform that I could do easily via pure SQL but just can't seem to get my head around it in the Field Calculator.

I have a table of 5000 weeds of which there are around 10 types.
I want to update one field in the table (called CLASS) to be either A, B, or C depending on another field: WEED_NAME

So, if WEED_NAME = 'Hyptis' then Class = 'A'
    else if WEED_NAME = 'Parkinsonia' then Class = 'B'  etc


texts: Best data provider for calculations in large vector files

texts:
I am building a QGIS plugin to read text files, convert them to vector files (100,000 to 1,000,000 features; ~5 numeric attributes) and display an additional attribute on a subset of the base vector file.

At the moment, I am loading the base QgsVectorLayer in the memory, then making a copy of a subset in the memory and adding the field/column.  The attribute to be appended is loaded into a numpy array, then updated in the copy.  Everything is done in the memory.  The process is slow though and I can see requesting and setting attributes are the main bottlenecks.

The additional data to be displayed is updated frequently (by the plugin itself).

Would it be faster to process and update the data on an external data provider (e.g. text file with WKT, or PostGIS), setting the attribute there, and then refreshing the source in QGIS?  Or is the QGIS memory faster than external processing + source refreshing?

(Note: the attribute requests and updates have already been optimized to operate on the least possible amount of data).


texts: How to merge (transfer) one layer into another PostGIS layer?

texts:
I want to merge two layers which have identical structure.

First layer is from postgres database table say "layer1" and  other layer is shapefile say "layer2". I want to merge layer2 "INTO" layer1 (NOT third layer).

I tried to copy paste in editing mode but it does not add attributes. any suggestions?

Note: the main layer is layer1 -postgresql database table.


texts: Accessing secure Geoserver layer with username password

texts:
I've created Users, Groups and Roles in Geoserver. And I set the user roles for my workspace.

user1 can read
user2 can write

And I want to access my layers as anonymous like this:

http:localhost:8080/geoserver/wms?request=GetCapabilities

This request returns 404 not found page.

I am using username password parameters but not working.

http:localhost:8080/geoserver/wms?request=GetCapabilities&user=some&password=some

If I login Geoserver admin panel and layer preview, I am accessing layer. If I logout, can not.


texts: Convert list of number to list of strings in Google Earth Engine

texts:
Is it possible to convert a list of numbers to a list of string in Google Earth Engine? I'd like to take something that looks like this

var lis = ee.List([27141, 20043]);

and convert it to something that looks like this

var lis = ee.List(['27141', '20043']);

I imagine this is really simple but I can't find a solution in GEE documentation.


texts: How to simulate/trigger a click on the map/view

texts:
In the ArcGIS Server JavaScript API v4x, is it possible to simulate or programatically trigger a click on the map at a given location?
I've seen Simulate a mouse-click inside an ArcGIS Online JavaScript Application as a workaround, but in my case I'd like to trigger the view's hitTest function, which runs after the user clicks. This is a shorthand way to get all visible features on the map, without having to iterate through them all and manually run a QueryTask.
Edit: to illustrate why this is useful, I have a map with 7 layers. When I manually click on the map, the view's hitTest function instantly returns the attribute values of all features at the location clicked. See the response.results object in the debugger:

There was no need for me to run 7 asynchronous QueryTask operations then collate the responses after waiting for each to return - the view handles this without requiring a single network request. The view also handles changing layer visibility automatically, so if (eg) one layer was switched off when I clicked, then only 6 graphics would be returned, etc.
Given a known coordinate, is it possible to trigger a click at that location?


texts: Problem with localhost:8080 using WMS and WFS services

texts:
I am facing a problem and I believe that my localhost port may be the main issue.

First of all, I am trying to connect my QGIS with these links of WMS
 (http://localhost:8080/geoserver/ows?service=WMS&version=1.3.0&request=GetCapabilities&) and WFS (http://localhost:8080/geoserver/ows?service=WFS&acceptversions=2.0.0&request=GetCapabilities&) services.

However, at my PC the connection fails, referring me that the connection is refused.

I think that the GeoServer service which I have installed two day ago, might be the main problem. I connected it with localhost:8080, reserving this port close. As a result, the links of the services can't connect.

I tried to uninstall GeoServer, reinstalling it using a different port (1234), but the problem wasn't resolved. I think that the 8080 port is still reserved from the first installation.

Also, I have requested the services in another PC, using QGIS, and the services worked fine.
Does anyone know how to resolve this problem?


texts: Error Accessing Dockwidget ComboBox. QGIS Plugin. (pyqgis)

texts:
I'm writing a QGIS plugin but I'm facing some basic problems. The plugin is a dockwidget. I've built it with plugin Builder and then compiled the resources file.

When I try to access the comboBox to populate it with values, an error happens:

  self.dockwidget.comboBox.clear()          AttributeError: 'NoneType' object
  has no attribute 'comboBox'

I checked in the Qt _Designer the object name and it's comboBox. When I do the same with a dialog plugin, instead a dockwidget, everything seems to be right, and the comboBox is populated. This is the code for the run method.

def run(self):
    """Run method that loads and starts the plugin"""
    self.dockwidget.comboBox.clear()
    for i in range(9):
        self.dockwidget.comboBox.addItem(str(i))

    if not self.pluginIsActive:
        self.pluginIsActive = True

        #print "** STARTING PruebaExp"

        # dockwidget may not exist if:
        #    first run of plugin
        #    removed on close (see self.onClosePlugin method)
        if self.dockwidget == None:
            # Create the dockwidget (after translation) and keep reference
            self.dockwidget = PruebaExpDockWidget()

        # connect to provide cleanup on closing of dockwidget
        self.dockwidget.closingPlugin.connect(self.onClosePlugin)

        # show the dockwidget
        # TODO: fix to allow choice of dock location
        self.iface.addDockWidget(Qt.LeftDockWidgetArea, self.dockwidget)
        self.dockwidget.show()

I don't know if it's related, but when I try to create a signal this error is show:

  File "XXXXXXXXXXX/PruebaExp.py", line 176, in initGui
      self.dockwidget.signalBox.connect(self.enviadato) AttributeError: 'NoneType' object has no attribute 'signalBox'

Same error when I try to acces a pushButton:

  File "/xxxx/xxxx/.qgis2/python/plugins/PruebaExp/PruebaExp.py", line
  221, in run self.dockwidget.pushButton.clicked.connect(self.enviadato)
  AttributeError: 'NoneType' object has no attribute 'pushButton'

Maybe is due to something I missed during the building of the plugin. I'm using QGIS 2.18.14.


texts: Synchronize from QField (QFieldSync) Not Working - Getting Handle Bad Layers pop-up

texts:
I am using QField with QGIS to edit data in the field. I used QFieldSync to package my QGIS project and was able to open it on my Android device and edit/add features just fine. I then copy the folder from my Android device to my desktop; go to my QGIS project; and select the folder through QFieldSync's 'Synchronize from QField'. The result is the 'Handle Bad Layers' pop-up in QGIS and I don't know how to proceed.

Can anyone advise on what I'm doing wrong? Here is some information:

Data: .shp files

QFieldSync Project Configuration: I've set the Action for all three layers to offline editing since I need to be able to make changes to the .shp files in the office while others are making updates in the field.

Python version: 3.6.0

QGIS version: 3.0.3-Girona

QField version: QField for QGIS - RELEASE CANDIDATE 3


texts:
I want to create a HillShade from a STRM dataset (spatial resolution: 30m) and I realized it has some noise on it (Little squares in the whole image).
I'd like to use Smooth Filters. I presume that QGIS has filters implemented. Could someone explain me how to use it and where I find it?


texts: How to access the datasource from a REST service

texts:
I have a ArcGIS REST service, is it possible to access the data source of the layer that is within the rest service. For example I have a layer called "Conservation Area" within the service, that references the dataset "PLANNING.CONSERVATION_AREAS". It is "PLANNING.CONSERVATION_AREAS" that i want to retrieve back. Is this possible

EDIT
I want to get access to the actual feature class name which is held within the database (oracle in my case).
I have read about the Data Source object REST Data Source Object which does suggest it is possible ? (but not sure how).

Thanks


texts: How to fix pgr_dijkstra error, query must return columns 'id', 'source', 'target' and 'cost'?

texts:
I followed the steps in "pgRouting 2.0 for Windows quick guide". At the step of running pgr_dijkstra using the sql below:

SELECT id, id1 AS node, id2 AS edge, cost, the_geom
  FROM pgr_dijkstra(
    'SELECT gid, source, target, st_length(geom) as cost FROM public.roads_svy21',
    1, 3000, false, false
  ) as di
  JOIN public.roads_svy21_vertices_pgr pt
  ON di.id2 = pt.id ;

It returned me an error message below:

Error, query must return columns 'id', 'source', 'target' and 'cost'

What is the possible cause of the error?


texts: Moving dead ahead: rhumb line, great circle or none

texts:
I think I understand the whole rhumb line and great circle idea but I can't seem to find an answer on this question. Suppose I am driving a car (i don't use ship because I don't want to put the constant compass bearing idea on the table) and that I can move for quite a long distance on land (across Russia for example) without any obstacles. If I keep the steer fixed so that the car moves dead ahead will I be following a rhumb line, a great circle or none of these.

I've found this on wikipedia:

  If you were to drive a car along a great circle you would hold the steering wheel fixed, but to follow a rhumb line you would have to turn the wheel, turning it more sharply as the poles are approached

but everything else I've found seems to contradict this. And yet I'm not sure if keeping the steer fixed means fixed dead ahead or fixed turned at 30 degrees left for example.

EDIT
After a request to give some examples of contradicting information found on the internet I edit my question to quote an example:

http://www.maritimeprofessional.com/Blogs/Maritime-Musings/February-2011/Rhumb-line.aspx

  The problem for navigation is that sailing a true great circle would
  require constant changes in the course steered

I now understand that what the writer of this text probably refers to when writing course steered is holding the compass indication unchanged. But this way of phrasing it is (at least to a not native english speaker like me) very confusing. In my mind "sailing a true great circle would require constant changes in the course steered" means that i have to constantly turn the ship's steering wheel in order to reach my destination.

I had found some similar examples which I can't find again now. But the greatest problem was I couldn't find any other source stating clearly that moving dead ahead equals following a great circle.

However, stating that "everything else I've found seems to contradict this" was a bit of an exaggeration and I apologize for it.


texts:
i implemented a compass for my three.js scene based on https://jsfiddle.net/prisoner849/gpbyw3co/ , but when i rotate object in my scene,the compass is is not showing the north correctly. Now i search for the cause. Could it be:
1. that i use TrackballControls
2. on camera lookAt() is called
this.camera.position.set(7000, 3000, -1750);
this.camera.lookAt(new Vector3(0, 0, -1750));

3. on TrackballControls target.set() is called
this.controls.target.set(0, 0, -1750);

this.controls.rotateSpeed = 4;
this.controls.staticMoving = true;
this.controls.noZoom = true;

4. zoom to mouse pointer is implemented
zoomToMousePointer(event) {
    event.preventDefault();
    event.stopPropagation();

    var vector = new Vector3();
    var dir = new Vector3();
    var position = this.camera.position;

    var rect = this.renderer.domElement.getBoundingClientRect();
    var left = event.clientX - rect.left;
    var top = event.clientY - rect.top;

    this.camera.updateMatrixWorld(true);
    vector.set((left / rect.width) * 2 - 1, -(top / rect.height) * 2 + 1, 0.5);
    vector.unproject(this.camera);

    var dir = vector.sub(position).normalize();
    let scalar = event.wheelDelta * 4;

    var vector = new Vector3(dir.x * scalar, dir.y * scalar, dir.z * scalar);
    this.camera.position.addVectors(this.camera.position, vector);
    this.controls.target.addVectors(
      this.controls.target,
      vector
    );
    this.renderer.render(this.scene, this.camera);
  }

6. Camera up manipulated
camera.up.set(0, 0, 1);

I use Three.js version 0.105.2
Seems to be the TrackballControls, when i change the fiddle to TrackballControls, compass is not working well anymore:


texts:
filter = 'LinkID IS NOT NULL'
cur1 = arcpy.UpdateCursor("ODDS_Links", filter)

high_id = 0
for row1 in cur1:
 if high_id < row1.LinkID:
  high_id = row1.LinkID

filter2 = 'LinkID IS NULL'
cur2 = arcpy.UpdateCursor("ODDS_Links", filter2)

i = high_id
for row2 in cur2:
 i+=1
 row2.LinkID = str(i)
 cur2.updateRow(row2)

The error I am getting is below. I know this works without the counter but something is wrong with counter and it won't add by one. If I comment out the counter it displays the error. I know I should use da but really don't feel like changing things around when I am so close.

  TypeError: coercing to Unicode: need string or buffer, int found


texts: Calculating size of area where data values exceed certain limit

texts:
I am using Python and Basemap to plot sea ice concentration within the Arctic region. I need to calculate the total area (square meters) where the sea ice concentration exceeds a certain limit (e.g. 15%), but I am unsure as how to do that.

My code for plotting the map is:

levels=np.arange(np.min(mydata),np.max(mydata),5)
mymap = Basemap(projection='npstere',lon_0=0,boundinglat=65)
x, y = mymap(lon,lat)

CS1 = mymap.contourf(x,y,mydata,levels,
                   cmap=mpl_util.LevelColormap(levels,cmap=cm.RdBu_r),
                   extend='max',
                   antialiased=False)

mymap.drawparallels(np.arange(-90.,120.,15.),labels=[1,0,0,0]) # draw parallels
mymap.drawmeridians(np.arange(0.,420.,30.),labels=[0,1,0,1]) # draw meridians
mymap.drawcoastlines()
mymap.drawcountries()
mymap.fillcontinents(color='grey')

My guess is that I need to convert the contours to polygons and then calculate the areas.


texts: Exporting query to identify what polygon certain points fall in using ArcGIS Desktop?

texts:
I'm using ArcGIS Desktop 10.

I have shapefiles consisting of points (street addresses).

I also have shapefiles consisting of polygons (ie, congressional districts, city council districts, etc... a separate shapefile for each kind of district)

How do I run a query that exports a csv of the address and what district that address falls in?


texts: How do I assign a height to each individual contour line using QGIS. Raster > Extraction > Contour

texts:
What I'm essentially asking is, is there a way to, when making contour lines using extraction in QGIS, assign heights to each line? For example instead of 1000, 2000, 4000, 6000 etc. etc. you have 500, 1000, 2000, 4000, 7000, 10000. So I want to put the lines at pre-chosen heights.
For example in this map it doesn't go up by the same amount each time, it goes up by a little at first and then slowly more. How do I do that in QGIS?


texts: Convert Arcpy results to feature class in ArcMap 10.1

texts:
I am using ArcGIS 10.1 and am a newbie in Python. Currently we are executing stored procedure (from sql server) through ArcPy in ArcMap. I want to create feature class from executing following python script.

I researched on Google but unable to find proper way/parameter in ArcPy. Currently going through following parameters GetParameter & GetMessage

Any help/steps will be great!

My ArcPy Script :

import sys
import arcpy
from arcpy import env

try:
    # Make data path relative
    arcpy.env.workspace = sys.path[0]

# Two ways to create the object, which also creates the
# connection to ArcSDE.
# Using the first method, pass a set of strings containing
# the connection properties:
#   <serverName>, <portNumber>, <version>, <userName>, <password>
#   arcpy.ArcSDESQLExecute("locolhost","5151","#","Admin","password")
# Using the second method pass the path to a valid ArcSDE connection file
sde_conn = arcpy.ArcSDESQLExecute(r"C:\Users\Administrator\AppData\Roaming\ESRI\Desktop10.1\ArcCatalog\test.sde")

# Get the SQL statements, separated by ; from a text string.
#sql_statement = 'select Name from all_points'
sql_statement = 'Exec [dbo].[sp_test]'
sql_statement_list = sql_statement.split(";")

print("+++++++++++++++++++++++++++++++++++++++++++++\n")

# For each SQL statement passed in, execute it.
for sql in sql_statement_list:
    print("Execute SQL Statement: {0}".format(sql))
    try:
        # Pass the SQL statement to the database.
        sde_return = sde_conn.execute(sql)
    except Exception as err:
        print(err)
        sde_return = False

    # If the return value is a list (a list of lists), display
    # each list as a row from the table being queried.
    if isinstance(sde_return, list):
        print("Number of rows returned by query: {0} rows".format(
            len(sde_return)))
        for row in sde_return:
            print(row)
        print("+++++++++++++++++++++++++++++++++++++++++++++\n")
    else:
        # If the return value was not a list, the statement was
        # most likely a DDL statment. Check its status.
        if sde_return == True:
            print("SQL statement: {0} ran successfully.".format(sql))
        else:
            print("SQL statement: {0} FAILED.".format(sql))
        print("+++++++++++++++++++++++++++++++++++++++++++++\n")


texts:
Is it possible to create a new service definition file from an existing ArcGIS service?

I have a python script that updates data on a daily basis.  I have figured out how to restart/overwrite my existing service so that it will represent the change in the data.  However doing so uses the default settings for the new service.

If possible, I would like to create a new service definition file from an existing service so that I do not have to go back manually and change the following parameters and capabilities:

Max # of Records Returned by Server,
Lock Database,
Max # of instances per machine,

Any advice would be great.  Thanks

Here is the solution:

    # Deserialize response into Python object
    dataObj = json.loads(data)
    httpConn.close()

    # Edit desired properties of the service
    dataObj["minInstancesPerNode"] = 1
    dataObj["maxInstancesPerNode"] = 10

    newdict = dataObj["properties"]
    newdict["maxRecordCount"] = 5000

    # Serialize back into JSON
    updatedSvcJson = json.dumps(dataObj)


texts: Least Cost Paths or routing: does anyone have examples of novel costs?

texts:
By novel costs I mean not simply slope or tolls but perhaps more intangible values like safety, social conformity or aesthetic appeal?
My interest is in trying to understand why certain routes were preferred by different groups of people in nineteenth century England.
An obvious modern analogy would be choosing an evening route home to avoid particular rough pubs and bars.


texts: Adding table of certain attributes to map layout?

texts:
I am using ArcMap 10.2 and I'm looking for a way to add  a table containing certain fields from select elements from a featureclass to my map. Ive tried the add table option from the attribute table, but it doesn't quite give what I want (the table must be editable to look a certain way).

Are there any options out there for 10.2 like arcpy.mapping.tableview or will upgrading to 10.3 so I can get access to the "graphic table element" function be best?


texts: Finding the number of points that lie within a radius of another point

texts:
I'm trying to find the number of restaurants that lie within 500m of a grocery store. I'm so confused as to how to do this. Using the near tool, it seems to only find the closest one. I tried generating a near table, but I was having trouble understanding that the different columns were.

Sometimes there are grocery stores closer than 500m apart, so I only want to count the restaurants that are closest to that grocery store.


texts: How to take 2 points from One featureClass, create a line between those points, and add the line to a new featureClass (aka Shape File)?

texts:
I'm a total noob/beginner in ArcObjects. I'm coding in C#.
I want to extract 2 points from a given featureClass, create a straight line between those 2 points, create a new feature class and add those points to the new featureClass.
The following is code (I have taken a lot of code from internet), but I am facing problems.

some_function()
{
        IWorkspaceFactory workspaceFactory = new FileGDBWorkspaceFactory();
        IWorkspace workspace = workspaceFactory.OpenFromFile(@"myGDB1.gdb", 0);

        IFeatureWorkspace featureWorkspace = workspace as IFeatureWorkspace;
        IFeatureClass fc_tenActors = featureWorkspace.OpenFeatureClass("ten_points");

        ISpatialReference lineFeatureSpatialReference = GetSpatialReference(fc_tenActors);
        IFeatureClass fc_line1 = IFeatureWorkspace__CreateFeatureClass(featureWorkspace, "Line1", lineFeatureSpatialReference);
        fc_line1 = featureWorkspace.OpenFeatureClass("Line1");
        IFeatureCursor featureCursorInsert = fc_line1.Insert(true);
        IFeatureBuffer featureBufferInsert = fc_line1.CreateFeatureBuffer();

        IFeature pFeature1 = fc_tenActors.GetFeature(3);
        IGeometry pFirstPointShape = pFeature1.Shape;
        IPoint pFirstPoint = pFirstPointShape as IPoint;

        IFeature pFeature2 = fc_tenActors.GetFeature(4);
        IGeometry pSecondPointShape = pFeature2.Shape;
        IPoint pSecondPoint = pSecondPointShape as IPoint;

  //************* PROBLEM HERE************************************
        IFeature pLine = CreateSegment(pFirstPoint, pSecondPoint);
 //************* PROBLEM HERE**************************************
  //  This is problem. pLine is empty after this statement.

        featureBufferInsert.Shape = pLine.Shape;

        AddFields(featureBufferInsert, pLine);
        featureCursorInsert.InsertFeature(featureBufferInsert);
        featureCursorInsert.Flush();

}

public static IFeature CreateSegment(IPoint pFrom, IPoint pTo)
    {
        Line line = new LineClass();
        line.PutCoords(pFrom,pTo);
        return line as IFeature;
    }

private static void AddFields(IFeatureBuffer featureBuffer, IFeature feature)
    {
        // Copy the attributes of the orig feature the new feature
        IRowBuffer rowBuffer = (IRowBuffer)featureBuffer;
        IFields fieldsNew = rowBuffer.Fields;

        IFields fields = feature.Fields;
        for (int i = 0; i <= fields.FieldCount - 1; i++)
        {
            IField field = fields.get_Field(i);
            if ((field.Type != esriFieldType.esriFieldTypeGeometry) &&
                (field.Type != esriFieldType.esriFieldTypeOID))
            {
                int intFieldIndex = fieldsNew.FindField(field.Name);
                if (intFieldIndex != -1)
                {
                    featureBuffer.set_Value(intFieldIndex, feature.get_Value(i));
                }
            }
        }
    }


texts: Iterative Interpolation w/ Built-In Python Console Not Changing Columns in PyQGIS?

texts:
I am trying to iteratively interpolate a shapefile using the QGIS IDW interpolator.  What I need to do is take a point shapefile with time series data in the columns associated with each point, interpolate that column of data to raster, then move to the next column, creating a new raster for each column.

The process I've generally been following is outlined here:
How to compute an interpolation raster from the python console in QGIS?

I've modified a few things and have been working on a test case (since the input data hasn't been provided yet - format should be similar except ~300 iterations), but I can't get the interpolator to move past the first, 0-indexed, column.

Code runs without error, outputs the 4 test raster interpolation .asc files, and adds them to the QGIS layers panel, but they are all the same.  You can see in the layers panel snippet below that the raster cell bounds are all the same and visually if I toggle the layers on/off, they are all the same.

Code below (only change is that I modified actual folder paths to dummy variables):

import qgis.core
import qgis.analysis
import os

pathToFile = r"input_folder\interp_test.shp"

interp_layer = QgsVectorLayer(pathToFile, 'interp_test','ogr')

#not used yet
field_names = [field.name() for field in interp_layer.pendingFields()]

#test first 4 columns
for i in np.arange(0, 4):
    layer_data = qgis.analysis.QgsInterpolator.LayerData()
    layer_data.vectorLayer = interp_layer
    layer_data.zCoordInterpolation=False

    #this line should change the interpolation attribute/column, but it doesn't
    layer_data.InterpolationAttribute = i
    layer_data.mInputType = 1 #points

    #pick interpolator engine
    idw_interpolator = qgis.analysis.QgsIDWInterpolator([layer_data])

    export_path =r"output_folder\rasterout_{}.asc".format(i)
    rect = interp_layer.extent()
    res = 500
    ncol = int( ( rect.xMaximum() - rect.xMinimum() ) / res )
    nrows = int( (rect.yMaximum() - rect.yMinimum() ) / res)

    output = qgis.analysis.QgsGridFileWriter(idw_interpolator,export_path,rect,ncol, nrows,res,res)
    output.writeFile(True)
    iface.addRasterLayer(export_path, "interpolation_output_{}.asc".format(i))

Input test attribute table:

Output to QGIS console:

Running 64-bit QGIS 2.18.3 on windows, using the built-in python console.


texts: cesium.camera.flyTo can't be called twice inside a  tag

texts: Is cesium.camera.flyTo() function a deferred javascript function? I mean this function is executed only after <script> is parsed. So what if I have to call cesium.camera.flyTo() more than once inside my <script> tag? The second call overrides the first call. Is there any better to get flights to multiple destinations in one go - like a loop of flyTo's or something?


texts: Localized Linear Histogram Match using python

texts:
I would like to use the local linear histogram match for filling the missing pixels of a Landsat 7 image with another older image (SLC Gap-Filled Products).

Do you have any idea how can I do it by using the gdal commands?


texts: Is there a way to set the data driven scale from arcpy? Looking at the documentation for the DataDrivenPages class, it doesn't seem to be accessible from outside of ArcMap. The script I run within the current map document exports the maps as jpgs based on the Name and Page Number fields I set, and at the Data Driven Scale I specified. When I put the same script into a model (with the necessary adjustments made) and run it from ArcCatalog, I get the same result, except the maps are now at Best Fit scale, which is the default for data driven pages.


texts: Deactivating editing toolbar for multiple features

texts:
I am creating a toolbar for a number of features as follows:

 for (var graphicNum = 0; graphicNum < graphicList.length; graphicNum++) {
        if (graphicList[graphicNum].symbol != undefined) {
            //Initializing the toolbar to be used for the feature modification process
            moveToolbarLabel = new esri.toolbars.Edit(app.map);
            moveToolbarLabel.activate(esri.toolbars.Edit.MOVE | esri.toolbars.Edit.SCALE | esri.toolbars.Edit.ROTATE, graphicList[graphicNum]);
        }
    }

The features all get the scale and rotation points around them as expected. However when I click a button to remove these scale and rotation points by deactivating the toolbar as follows:

if (moveToolbarLabel != undefined) {
        moveToolbarLabel.deactivate();
    }

the scaling and rotation points do no go. Is there any additional code that has to be included since there would be multiple toolbar due to having the toolbar activated for multiple features?


texts: python cluster Points with DBSCAN keeping into account Polygon boundaries

texts:
I have a GeoPandas dataframe called geo containing:

431978 unique spatial Points representing all households in a city.
112737 Polygons representing all the premises (panden) in a city

dataframe looks like this (note that the x, y coordinates come from the Point geometry column):

I need to cluster all these Points within a 5 meter range using DBSCAN but partitioned over the premises so that all clustered Points cannot cross the boundaries of the Polygon.

So the following code produces the WRONG clustering as theoretically depicted in the screenshot:

from sklearn.cluster import DBSCAN
epsilon = 5
coords = geo[['x', 'y']].values
db = (DBSCAN(eps=epsilon, min_samples=1, algorithm='ball_tree', metric = 'euclidean')
.fit(coords))

cluster_labels = db.labels_
unique_labels = set(cluster_labels)
# number of clusters
num_clusters = len(set(cluster_labels))
print ('number of clustered verblijfsobjecten:
{}'.format(num_clusters))

The right DBSCAN cluster assignment, again, would take into a count the boundaries of the Polygons that the clustered Points CANNOT cross. So to depict the right clustering for one of the premises:

I assume assigning the DBSCAN algorithm on each group resulting from the geo.groupby(['pand_id']) and get_group method, would be a way to go ?

In SQL this is easy but my question is how to do this in Python/sklearn. thnx!

SQL statement:

Q = """drop table if exists TABLE.SCHEMA;
create table service_afvalcontainers.stag_vo_cluster as
select
    pand_id,
    ST_X(geometrie) x,
    ST_Y(geometrie) y,
    ST_ClusterDBSCAN(geometrie,5,1) over(PARTITION BY pand_id) as
    cluster,
    geometrie
from
    TABLE.SCHEMA bla bla

;
"""
UPDATED ANSWER PART: In Python I tried:

### get pand_id to feed into loop.
unique_pand_id = geo.pand_id.unique().tolist()

coords_frame = []
for pand in unique_pand_id[:10]:
coords = (geo.loc[geo['pand_id'] == pand, ('index' ,'pand_id' ,'verblijfsobject_id', 'x', 'y')].values)
coords_frame.append(coords)

cluster_label_frame = []
epsilon=5
for elem in [[elem[3:5] for elem in group] for group in coords_frame]:
db = (DBSCAN(eps=epsilon, min_samples=1, algorithm='ball_tree', metric = 'euclidean')
 .fit(elem))
cluster_labels = db.labels_
cluster_label_frame.append(cluster_labels)

This cluster_label_frame looks ok indeed but the thing is, I need to have these cluster_ids and attach them back onto the original geo dataframe. I probably need to add the geo index to this label_frame with which I can hopefully join back to the geo frame. It doens' t feel like the right approach but it`s a step further...


texts: Direct connecting to PostgreSQL from ArcGIS Desktop?

texts:
I am trying to connect to Postgresql database that I just installed on my server machine. As instructed by Esri I have download from postgres and installed psqlodbc. However, when I tried to add the database connection using "Add OLE DBconnection" I could not find it listed under OLE DB provider(s). More over I do not understand the connection parameters. I could the see driver installed in the control panel (under add/delete programs). However I could not see it listed in the OLE DB connection. Can you help me understand my problem? Thanks.

Additional information: I just install postgresql. It has a default database. I have not install ArcSDE on top of it.


texts: Empty shapefile created with Python and GDAL

texts:
I need to create a point shapefile using the coordinates specified in the variable I called "my_points_dataset".

'my_points_dataset' variable has this structure:

1615311.055743243079633 4820229.073873873800039 0,1615321.055743243079633 4820229.073873873800039 0,1615331.055743243079633 4820229.073873873800039 0,1615341.055743243079633 4820229.073873873800039 0,1615351.055743243079633 4820229.073873873800039 0, etc`

The class of this variable is 'osgeo.ogr.Geometry'

Each point in 'my_points_dataset' is defined as x, y and z coordinates (z could be negligeable because is always 0).
x, y and z data are separated by a space, each point is separated by comma.

I wrote this Python code, but at the end the result is an empty shapefile.
Do yo know why?

import ogr, gdal
shpDriver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists('/path/mypoints.shp'):
        shpDriver.DeleteDataSource('/path/points_shapefile.shp')
    outDataSource = shpDriver.CreateDataSource('/path/points_shapefile.shp')
    outLayer = outDataSource.CreateLayer('/path/points_shapefile.shp', geom_type=ogr.wkbMultiPoint)
    featureDefn = outLayer.GetLayerDefn()
    outFeature = ogr.Feature(featureDefn)
    outFeature.SetGeometry(my_points_dataset)
    outLayer.CreateFeature(outFeature)
    outFeature = None


texts:
I have location data points for colonies (csv file). I have a shapefile that I need to duplicate and georeference to each colony. Is there a way to do this in r?

EDIT: I want to create a copy of the existing KUD shapefile for each colony location. I do not want to create a unique KUD for each colony.

Picture shows the shapefile (green) that I am trying to replicate for each colony and have the center of the shapefile georeferenced to on the colony. Arrows show the change I want to make in the location of the shapefile.

Picture shows the data table of the colonies. They are imported as a csv file and I need to create a copy of the shapefile for each colony data point located at the same lat/lon.


texts: PyQGIS QgsSvgMarkerSymbolLayerV2 setDataDefinedProperty() for color doesnt work

texts:
I have a script which works well and changes the symbol QuarterCircle's color and angle :
vLayer = QgsVectorLayer(uri, "myData" , "delimitedtext")

symbol = QgsSymbolV2.defaultSymbol(vLayer.geometryType())

simple_marker = QgsSimpleMarkerSymbolLayerV2()
simple_marker.setShape(QgsSimpleMarkerSymbolLayerBase.QuarterCircle )
simple_marker.setSize(6.0)

ddpAngle=QgsDataDefined( True, True, "Direction")
ddpColor= QgsDataDefined( True, True, "CASE WHEN DATA % 3 = 0 then  color_rgb( 255,0,0) WHEN DATA % 3 = 1 then  color_rgb( 0,255,0) WHEN DATA % 3 = 2 then  color_rgb( 0,0,255) END" )
simple_marker .setDataDefinedProperty( "angle", ddpAngle)
simple_marker .setDataDefinedProperty( "color", ddpColor)

symbol.changeSymbolLayer(0,simple_marker)

vLayer.rendererV2().setSymbol(symbol)
vLayer.triggerRepaint()

I wanted to change the QgsSimpleMarkerSymbolLayerBase.QuarterCircle symbol to user defined .svg format symbol.
symbol = QgsSymbolV2.defaultSymbol(vLayer.geometryType())
svg_marker = QgsSvgMarkerSymbolLayerV2( 'D:\\sample.svg' )
svg_marker.setSize(12.0)
svg_marker.setOutlineWidth(0.25)
ddpAngle=QgsDataDefined( True, True, "Direction")
ddpColor= QgsDataDefined( True, True, "CASE WHEN DATA % 3 = 0 then  color_rgb( 255,0,0) WHEN DATA % 3 = 1 then  color_rgb( 0,255,0) WHEN DATA % 3 = 2 then  color_rgb( 0,0,255) END" )
svg_marker.setDataDefinedProperty( "angle", ddpAngle)
svg_marker.setDataDefinedProperty( "color", ddpColor)

symbol.changeSymbolLayer(0,svg_marker)

vLayer.rendererV2().setSymbol(symbol)
vLayer.triggerRepaint()

For this one it changes the angle  but the color part doesn't work and all the colors are black.
Can you please help me find the problem?
UPDATE: I realised that when I do it by manually from desktop app it works(entering data defined override by manually). But when I paste my code to desktop app's python console by changing vLayer to iface.activeLayer() it also doesn't change color and it shows no data defined override for color part of it, it just seems deactivated.
Can anyone give an example portion of code that imports an svg and modify its color according to some data belongs to a vector layer?


texts: circle with dfferent colours representing different criteria

texts:
I want to draw circles on maps with different colours, this gives the signal strength of the tower, suppose circle is of 50 m then first 20m should be given a color1 and next 20 m color2 and remaining color3.

Is this possible in Geoserver SLD or openlayers? How?


texts: Error Message when Running Global Morans with Spatial Weights Matrix

texts:
I am attempting to run a global Moran's I test with a spatial weights matrix, however when I run the process I always receive an error message and I was hoping someone could point out when I am going wrong.

My project is looking at spatial autocorrelation of health in London, as part of the process was looking to produce a spatial weights matrix using Queen's method and K nearest neighbors in order to produce a global Moran's statistic using the matrix. However I always receive the error message attached.

How do I need to set up the spatial weights matrix to avoid these errors?


texts: PyQGIS python.exe has stopped working but my output is still produced

texts:
I am running a script in QGIS which loops around a set of rasters performing a number of different processes ultimately outputting a new raster and a vector per input. At the end of each loop i get the error message 'python.exe has stopped working' but my output is still created and saved in the destination folder.

The script I am using is as follows:

##Input_raster=multiple raster
##result_shp=folder

import os

layers = Input_raster.split(';')
for raster in layers:
    #Names for output
    result_name = os.path.splitext(os.path.basename(raster)) [0] + "_Extent.tif"
    shp_name = os.path.splitext(os.path.basename(raster)) [0] + "_Outline.shp"

    #Take grid and create extent
    out1=processing.runalg('gdalogr:rastercalculator', raster,'1',None,'1',None,'1',None,'1',None,'1',None,'1','A>0',None,5,None,result_shp + "/" + result_name)
    #Ratser extent to polygon
    out2=processing.runalg('gdalogr:polygonize', out1['OUTPUT'],'DN',None)
    #Calculate area m^2 and add to new field 'Area_m2'
    out3=processing.runalg('qgis:fieldcalculator', out2['OUTPUT'],'Area_m2',0,10.0,3.0,True,'$area',None)
    #Select polygons >5m^2 and extract to new shp
    out4=processing.runalg('qgis:extractbyattribute', out3['OUTPUT_LAYER'],'Area_m2',3,'5',None)
    #Fill holes less than 200m^2
    out5=processing.runalg('qgis:fillholes', out4['OUTPUT'],200.0,result_shp + "/" + shp_name)

Is there a reason this might be occurring for each loop?

Do i need to somehow 'close' the input or output at the end of each loop to avoid this happening?

I have now run the same process using QGIS Graphical Modeller and found the same error stopping the process at step 1 GDAL raster calculator. When I view the log the final message given is 'GDAL command output:'


texts: Get values from non spatial table in geopackage via Python to compalist

texts:
I need to get values from non spatial table in geopackage via Python. Do You have any ideas how to get this?
I'm trying get this using spatialite3 library and simple string select query but I have receive error message "there is no table_name in database"
When I'm trying get this from qgis layer list I can't find appropriate way to do this. I'm trying this:

database = r'path_to_my_geopackage.gpkg'
uri = QgsDataSourceUri()

uri.setDatabase(database)
schema = ''
table = 'tablename'
geomColumn = ''

content = uri.setDataSource(schema,table,geomColumn)


texts: Problems with spatial autocorrelation in Geographically Weighted Regression (GWR) analysis

texts:
I have conducted an Ordinary Least Squares (OLS) and Geographically Weighted Regression (GWR) analysis in ArcGIS where I have tried to predict biodiversity patterns in an area. When I examine the results with Spatial Autocorrelation the residuals tends to be clustered with a fairly high Z-value.

The result from the GWR, however, tends to be the expected when I visualize it. How do I handle the clustering? What problems does it creates, IS there any problems?


texts:
Can anyone lend a hand with what I'm doing wrong here? I'm trying to get a 92x92 (because it's the biggest I can seem to setup) PDF to come out of ArcMap 10.1. I'm in layout view, set 92x92 in the page size, set it to scale map elements proportionally, etc.

It comes out not actually filling the whole page.

Ideas?


texts: View Image instead of label in Menu on arcmap addin toolbar

texts:
I want to set image on menu in ArcMap addin toolbar. is there any way to set image ? I have tried but this is not working..
My code is here in config.esriaddinx:

--Menus>
        Menu id="GridViewMenu" caption="Element View"  image="Images\gridview" isRootMenu="false">
          Items>
            Button refID="Button1" />
          /Items>
        /Menu>
/Menus>


texts:
I am developing a stand alone python script that will invoke methods from processing module.
I am seeing error "ImportError: No module named processing" on ubuntu.

Here is the code.

import PyQt4.QtGui
import PyQt4.QtCore
import PyQt4.Qsci

from qgis.gui import *
from qgis.core import *
import os

import processing

print 'hello'


texts: Moving all merged layers at once using QGIS

texts:
I am using QGIS 3.10.
Following the query at QGIS 3.10 problem with merging layers & multipart promotion I merged my layers into one piece. Now I would like to make them draggable across the map, as per my previous query at Moving layer across map canvas.
I spotted, that despite the all-in-one performance, the draggable tool still moves only one layer selected:

I want to have all of them moved at the same time.
Is it possible?

I tried to use the " Remove duplicate geometries" option, although the result is still the same:

A similar question has been asked as Moving two layers at once in QGIS but is not applied to the merged layers.
My attribute table looks as follows:


texts:
I am researching potential tech solutions for personnel to use out in remote areas in the event of an emergency situation.

I have been told I need to

1) find satellite based solutions as the phone network is not reliable.

2) Find a device that facilitates user check-in at set intervals so that an alert is raised if they fail to check in and an escalation procedure can begin.

I am struggling to find a single device that fulfills both these criteria. Can anyone make me recommendation or at least explain why these two criteria seem to be mutually exclusive in most devices. I know the check-in exists for phone apps, why can't it work for satellite based devices?


texts: Shapefile projection for Matplotlib Basemap

texts:
Preface: First time working with GIS data.

I'm using python's Basemap library for overlaying data on maps.

The tutorial I'm using includes a shapefile for London (london_wards.shp).

I use the bounds method to get the minimum bounding rectangle.

s="/data/london_wards.shp"
shp = fiona.open(s)
shp.bounds

(-0.5103750689005356,
 51.28676016315085,
 0.3340155643740321,
 51.691874116909894)

This is the latitude, longitude (51.5072° N, 0.1275° W) for London.

I pass these bounds to a basemap instance.

Next, I've done the same for shapefiles obtained from SF OpenData.

However, the bounds are apparently not latitude, longitude.

s="/data/sfpd_districts.shp"
shp = fiona.open(s)
shp.bounds

(5979385.3163340045, 2085856.3243659574, 6024751.466714004, 2131290.9014959573)

So, I can't pass them to basemap.

Several shapefiles from SF OpenData have similar units returned from shp.bounds.

What are these values?

How are they converted to latitude, longitude for further use?


texts: Saving separate output files when iterating ArcGIS model?

texts:
I need to run a model in ModelBuilder 1000 times in order to obtain the average distance of random points in a wildlife species' home range to anthropogenic features. I need to find the average of the random points from EACH of the iterations. However, when running the model 1000 times, the output file is continuously being replaced, so you only end up with one output table, the very last output table.

How can I use model tools in order to force the model to save separate files for EACH iteration so that I can find the average distance of the random points to the features for EACH iteration?

I cannot seem to find ANYTHING online about this problem!


texts:
I am trying to find an sld for a Ship. I am parsing an AIS feed and would like to make the ship different colors based on class of ship and length as well as rotate it based on the COG (course). Currently I have  an Arrow that rotates but it using a well known class. I would really prefer something that actually resembles a ship based on a polygon.

Here is the style I am using currently:

<PointSymbolizer>
  <Graphic>
    <Mark>
      <WellKnownName>arrow</WellKnownName>
        <Fill>
          <CssParameter name="fill">#FF0000</CssParameter>
        </Fill>
        <Stroke>
          <CssParameter name="stroke">#000000</CssParameter>
          <CssParameter name="stroke-width">2</CssParameter>
        </Stroke>
      </Mark>
      <Size>12</Size>
      <Rotation>
        <ogc:Sub>
          <ogc:PropertyName>COG</ogc:PropertyName>
          <ogc:Literal>90</ogc:Literal>
        </ogc:Sub>
      </Rotation>
    </Graphic>
  </PointSymbolizer>


texts: Unknown red rectangle appears in map window

texts:
I would like to know where this red rectangle comes from? It's not a polygon or a line, it can't be activated or removed or anything else. It just appeared during the project work (adding shp-files, georeferencing rasters and digitalizising,...)

Maybe it's about extents?!


texts: Why are grid and scalebar measures 1mm too long in a 1:10000 print?

texts:
I used qgis 1.8 to compose a map from the WMS-service of the Belgian NGI.  I added a scalebar and grid to it.  When I measure it by hand, the distance between the longitudelines is 101mm instead of 100mm.  The distance between the lattitudelines is correct.

Even the scalebar is 1mm to long.

Has somebody encountered the same problem yet?


texts:
I'm looking for a similar functionality as this tool: Import Metadata in ArcGIS Pro. I am upgrading some of our older scripts and I can't seem to find a similar tool or even one that could do something somewhat similar.
Does ArcGIS Pro have any ArcPy options for this at 2.6 (which I am currently using)?


texts:
Does anyone have any recommendations for outputting a PDF of map from Manifold 8?

The following printer drivers have all failed to output anything:

CutePDF
PDF995
eDoc

Nitro PDF worked (after 3 hours), but produced a 30mb mess of a map.  Primo PDF is still going, and has been processing for hours...

My preference is for the output to remain as vectors to keep the quality of the map. Not fussed if it's a GeoPDF or straight PDF.


texts: How can I use arcgis server to help me find a street according to Latitude and longitude?

texts:
I am now in a project using a taxi GPS data(latitude and longitude ) to locate the road segment it is being on.I have the road map data. So, I wonder if I can use arcgis server to get the road segment information according to the GPS data?I wander if there is an java api for me to call easily? for Example, I have the latitude 16.34211,longitude 4.0102,I call the function getRoadSegment(16.34211,4.0102),the function then return the road segment (maybe the startpoint(x1,y1) and the endpoint(x2,y2)) this point is in.!
This picture is when I import mxc data file into ArcGIS.


texts: How to save .shp to .csv so that the Name is in the first "column" of the .csv?

texts: I am trying to convert a point shapefile to .csv with locations in QGIS. The problem I am having with the conversion is mainly that, in the resulting .csv file I would like to have the Name and then Long, Lat, however it constantly comes out as Long, Lat, Name. Is there any way I can save it so that the Name is in the first "column" of the .csv?


texts: How to create a sub region shapefile from MasterMaps on QGIS?

texts:
I have the OS Mastermaps on a PostGIS database and this is connected to QGIS. Mastermaps is visualised well on QGIS with all the layer styles.

A region is clipped from PostGIS and created a new table. When this is visualised on QGIS the layer styles are not visible as shown in the image.

I tried copying the layer style from the original layer to the new layer but still it does not show on the layer. However in the properties the new style can be seen.

Can you please help me to apply the original styles to this new layer and save as a shapefile?


texts:
I'm a Windows user getting an error when using GDAL from the command line (which I'd like to do for a couple reasons).

Attempting to use gdal_polygonize.py (in my case), as in:

gdal_polygonize.py test\test-gdal-tmp.tif -f "ESRI Shapefile" test\test.shp test

results in the error:

from osgeo import gdal
ImportError: No module named osgeo

whereas the line "from osgeo import gdal" works fine if I just load up Python (e.g. some standard IDE).

Prompting the code to print out sys.path both in my IDE & from the GDAL command line actually results in a different set of paths. The IDE gives me the necessary & expected folders like C:\\Python27\\lib\\site-packages. The command line GDAL call actually gives me ARCGIS-related folders not in the IDE's path, e.g. C:\\Python27\\ArcGIS10.4\\lib\\site-packages & without the necessary path to find osgeo (obviously).

I'm sure I've done something sinful in installing GDAL &/or have some installation incompatibilities, but for the life of me I can't figure out what. How can I get command line GDAL not to change the system path?

Other info:

OS = Windows 10
To install GDAL, I grabbed the files from here: http://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal (GDAL‑2.1.3‑cp27‑cp27m‑win32.whl) and used pip
Have ArcGIS on my system (clearly)
Made some other attempts to install GDAL (e.g. with Anaconda) before realizing this wouldn't allow them to be integrated with other code as necessary for my project
Ultimately trying to get this code to work, hence the command lines: https://github.com/nypl-spacetime/map-vectorizer


texts: Why do I get Error 00728: Field "" Does not exist (occurs running AddField_management & CalculateField_management in succession)

texts:
I have a Python Toolbox and want to use the zonal statistics as table function on user-provided polygon shapefiles/feature classes that may or may not have multiple parts.

In the code, I add a new field to the user-provided polygon and create the zones (code below is abbreviated with only pertinent parts):

polygon = parameters[0].value
arcpy.AddField_management(polygon, "Zone", "TEXT", "", "", 15)

EDIT 1: (code immediately below)

d=arcpy.Describe(polygon)
if d.dataType == 'FeatureLayer':
     polygon = d.featureClass.catalogPath

Then using field calculator, I create sequential zone numbers:

expression = "autoIncrement(!FID!)"
codeblock = """def autoIncrement(fid):
    fid = int(fid)
    fid += 1
    return "Polygon " + str(fid)"""
arcpy.CalculateField_management(polygon, "Zone", expression, "PYTHON", codeblock)

When I try the tool, the CalculateField_management function stalls and I receive "Error 000728: Field "Zone" does not exist within table". However, opening up the shapefile used to test the tool does reveal that a new field called "Zone", string datatype and length of 15 has been created.

When I refresh the Python toolbox and rerun the tool, this time it runs correctly populating the rows ("Polygon 1", "Polygon 2", ...) and the ZonalStatisticsAsTable function is able to run.

I found a previous post on this issue: Why am I getting an ERROR 000728: Field *fieldname* does not exist within table even though field exists?

Attempted Workarounds:

(1) Instead of using arcpy.CalculateField_management, I tried using an update cursor.

arcpy.AddField_management(polygon, "Zone", "TEXT", "", "", 15)
start = 1
cursor = arcpy.da.UpdateCursor(polygon, "Zone")
for row in cursor:
   row[0] = "Polygon " + str(start)
   start = start + 1
cursor.updateRow(row)
del row, cursor

arcpy.sa.ZonalStatisticsAsTable(polygon, "Zone", ...)   [code goes on]

Instead of the previous error, the code stalls on the ZonalStatisticsAsTable function instead and I receive "Error 001000: Field "Zone" does not exist". As before, checking the shapefile's attribute table reveals that the "Zone" field exists as it should after arcpy.AddField_management function.

(2) Use "FID" as zones to avoid having to use arcpy.AddField_management and arcpy.sa.ZonalStatisticsAsTable in succession.

Tool does now run without either errors, yet summary statistics are incorrect. Basically it looks like values of "MIN", "MAX", and "MEAN" have been shifted over, so this doesn't look like a viable option.

My thoughts:
Is there some kind of lock preventing the code from running top to bottom? Issue occurs after I add the new field and then want to do something to that field. It's weird that I need to rerun the code twice in order to not receive errors. Since I hope to have this tool available for other users, this isn't acceptable.


texts: How can I web map service published by Dojo in Flex app

texts:
We want to include a base map layer in hosted Flex app. We requested web map service from the vendor and I got a file (.htm) from client with some dojo code. I can understand the code but I am not sure how can I include this in my flex app configuration:

    // info used to create web tile layers
    var mapurl  = "http://bnvnbvnn/${0}/${2}/${1}.png"
    var lyr = new esri.layers.WebTileLayer(mapurl);
    // clear out the side bar
    //dojo.byId("leftPane").innerHTML = "";
    // create and add the layers
    //var layerDiv = dojo.create("div");
    map.addLayer(lyr);

  }

If I will replace variables on URL, I can see images but how will I use this for Flex configuration- where it comes like this:

 <layer url="http://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer"
visible="true" type="tiled" label="ArcGIS Streets"/>

I have never worked with Dojo toolkit, so I have no idea (or may be there is some communication gap)

Thanks


texts: Weird regions in lagoon water ripples pattern (Google Maps)

texts:
When looking at some water reservoirs or lagoons from Google Maps (satellite), I've sometimes found this kind of structures showing regions with a slight color/bright difference:

They seem to be an artifact, not a physical effect on the ripples due to the wind:

However, the borders shape not always seem due to a mixture of different tiles (not always right angles):

I understand that they might need to use different sets of pictures for different parts of the coastline, but this does not seem to require the weird shapes to make everything consistent, especially when considering some of the regions that are not in contact with the coastline but just with other water regions. What is the reason for these shapes? Is this related to avoiding the image of the sun reflected on the surface? (every region exhibit a brightness gradient which may point to this image).


texts:
I'm currently attempting to export a selected set of data from a feature class in a database to another database as a new feature class. I understand everything perfectly right up to the IExportOperation.ExportFeatureClass bits. The last argument calls for an IFeatureClassName object. This appears (to me) to not be a path or a file name, but an object that maintains all of that data. I have been looking at some other questions that discuss or mention the method to export a feature class, but I cannot wrap my head around it properly.

What is an IFeatureClassName, and how does one create it for a file that does not exist yet (the Geodatabase does exist)?


texts: GeoTiff image not aligned after using gdal_translate and gdalwarp

texts:
I don't see what I am doing wrong. The coordinate system was assigned to the GeoTiff using the following gdal_translate command:

(Usage: gdal_translate -of GTiff -a_srs EPSG#### -a_ullr upper left lower right source.tif target.tif)

gdal_translate -of GTiff -a_srs EPSG:29191 -a_ullr 567.4 8321.8 571.4 8318.8 c5100102_pdf_clipped2boundary.tif output01.tif

followed by

gdalwarp -t_srs EPSG:29191 output01.tif final01.tif

This last step did not work.

Did I miss something?

I have gotten this far:

Something is still missing, as Tiff does not align perfectly with the vectors. Reprojecting rasters is a new area for me.


texts: Use of locked magnification dialog displaying USGS US Base Topo Tiled layer

texts:
I am studying the USGS Base Topo Layer for my railroad interest. I chose 1:10,000 magnification, as that renders the raster tiles optimal for sharpness. Larger magstep results in provider not giving data (grey tiles). When panning the tiles, the magnification strangely changes to 1:9999, 1:9998 etc. with each panning command, despite activating the lock magnification checkbox.
Is this a QGIS bug? I would think that locking the magnification at 1:10,000 should be honored by the system.
I am using QGIS-3.4.11-Madeira.


texts: Problem with Pygeoprocessing since upgrading to version 0.4.4 from version 0.2.2

texts:
I need to calculate D-Infinity Flow Direction on my input raster (geotiff).  I use pygeoprocessing for this.  Since I uninstalled the version that I had (0.2.2) and installed a newer version of the python package (0.4.4), I've not been able to calculate flow direction using pygeoprocessing.routing.flow_direction_d_inf module. I get the error message "AttributeError: 'NoneType' object has no attribute 'GetProjection'" which I am assuming means that it wasn't able to load the input image to get it's attributes.

This is what I type in at the Python Idle GUI:

>>> import os
>>> import sys
>>> import pygeoprocessing
>>> import pygeoprocessing.routing
>>> from osgeo import gdal
>>> pygeoprocessing.routing.flow_direction_d_inf('C:\\Projects\\data\\work\\DemPrime.tif', 'C:\\Projects\\data\\work\\DIFD.tif')
  File "C:\Python27\ArcGISx6410.1\lib\site-packages\pygeoprocessing\routing\routing.py", line 331, in flow_direction_d_inf
    elevation_raster_path_band[0], target_flow_direction_path)
  File "src\pygeoprocessing\routing\routing_core.pyx", line 711, in pygeoprocessing.routing.routing_core.flow_direction_inf (src/pygeoprocessing/routing\routing_core.cpp:12951)
  File "C:\Python27\ArcGISx6410.1\lib\site-packages\pygeoprocessing\geoprocessing.py", line 979, in get_raster_info
    raster_properties['projection'] = raster.GetProjection()
AttributeError: 'NoneType' object has no attribute 'GetProjection'

Has anyone else run into problems when upgrading pygeoprocessing?


texts: Why doesn't QGIS accept my manually altered geospatial column?

texts:
I have a table with a geospatial column with srid 28992 but noticed that if I try to show it in QGIS (v1.9) on any OpenLayers map (through the OpenLayers plugin) it is off by a few hundred meters.

Then I converted it to epsg:900913 (remembering something that with an old plugin this projection was required) and then it is shown on the correct spot.

However when I look in the metadata of this layer in QGIS (right-click -> Properties -> Metadata) everything is blank and when I do Zoom to Layer Extent the map jumps to water near Africa (suggesting it doesn't understand the projection).

Somehow the metadata gets lost to QGIS during my manual conversion. The script I use for the manual conversion is:

alter table <table_name> drop constraint "enforce_srid_<column_name>";
update <table_name> set <column_name> = ST_Transform(<column_name>,900913);
alter table <table_name> add constraint "enforce_srid_<column_name>" CHECK (st_srid(<column_name>) = 900913);
update geometry_columns set srid=900913 where f_geometry_column = '<column_name>';

The result of this seems no different (except that now I have an extra column) than if I use:

SELECT AddGeometryColumn('<table_name>','<new_column_name>',900913, 'LINESTRING', 2);
UPDATE <table_name> SET <new_column_name> = ST_Transform(<column_name>,900913);

Does any of you know what god particle I am looking for?


texts:
Is it possible to do SQL filters in ArcGIS online similar to the SQL tab of CartoDB maps?

I want to do an expression on several table fields.

Below is a screenshot of the map control panel that I have.


texts: NLCD Detection Limit of Land Classification

texts:
I am working on a project with NLCD 2016 data that involves calculating the percentage of each land cover class around sites. I have this working, but need to impute zeros around each site for the classes not measured around each site.
Part of this imputation involves knowing the 'detection limit' of the classification. I have looked quite a bit, but all I have been able to dig up is that the Terrestrial & Ecosystems Data from USGS has "a 30 meter pixel cell and in most areas a minimum mapping unit of 0.4 ha (1 acre) this means that small patches of vegetation can be missed in the modeling process."
However, this 1. isn't the data I'm using and 2. even if it was I am still not sure if the detection limit in this case would be 30 m^2 or 1 acre.
Looking for a detection limit (the minimum amount of a land usage class required to make a classification) for NLCD 2016 data.
Update: I should add that I'm like 99% sure what I want is the MMU. Just need to find it for the NLCD. I'm just a little confused why the MMU would be less than the resolution...


texts: PostGIS subquery make linestring from points

texts:
I have thousands of points and I am trying to intersect them with US states, sort them by date and then create a linestring from those points grouped by date, and state name. I have accomplished this by performing two separate queries, but I want to combine the process into one larger query. I am fairly certain this requires a subquery. here is what I have so far

create table mypoints as
select us.name,ralphpoints.* from us,ralphpoints where st_intersects(us.geom,ralphpoints.geom) order by ralphpoints.datereal asc;

create table lines as
select name as state,datetosort as date,st_makeline(geom) as geom from mypoints group by name,date order by date;

this query produces the below picture which is accurate.

Now when my subquery:

create table lines as
select t.name as state,t.datetosort as date,st_makeline(t.geom) as geom from(
select us.name,ralphpoints.* from us,ralphpoints where st_intersects(us.geom,ralphpoints.geom) order by ralphpoints.datereal asc)t
group by state,date order by date

the result from the subquery is off. i tried to just follow the same logic as the separate two queries, i am not sure why the results are different.

*us refers to US shapefile of states
*there are two date fields from the original ralphpoints table. 1.datereal = date day/month/year with detailed time - used for subquery order by to 2.datetosort = just day/month/year used in the group by


texts: Select polygons one by one and export them using arcpy

texts:
I have a shapefile that contains 200+ metropolitan cities. In each city, there are some tract polygons. The dbf table looks like:

tract  city    xvar
1      0040    20
2      0040    12
3      0040    45
...    ...     ...
5      5602    22
6      5602    44
...    ...     ...

My goal is to select a city and save it as a separate shapefile the name of which is the same as the city code:

0040.shp

tract  city    xvar
1      0040    20
2      0040    12
3      0040    45
...    ...     ...

5602.shp

tract  city    xvar
5      5602    22
6      5602    44
...    ...     ...

How should I automate this process for all the cities using arcpy?


texts: Summing attributes from different layers

texts:
I'm working on a project in QGIS and I'm stuck. I have some shapefile describing the distribution of an animal,  these file have an attribute name 'PRESENCE' with value 1. There are point when only one file is present and other where there are more.
My goal is to create a new layer (I think it's should be a raster) rappresenting the sum of all the 'PRESENCE' attributes for every point.
I've tried using the function join attributes by location but it doesn't seem working.
The image is showing how my files look and the attributes table

Can you help me?


texts:
I am new to R and would like to generate Latin Hypercube samples (LHS) for three variables within the following limts:

1.25 ≤ x1 ≤ 5 ,
1.25 ≤ x2 ≤ 5
7 ≤ x1 ≤ 25

with total 3*10 = 30 samples.

What is the syntax to be used in Rstudio?


texts: QGIS Assign point data to a line layer considering line direction

texts:
I have a line layer representing a river bank. The river bank is divided in numerous separate line features representing zones, each separate line has many vertices. The lines have a set direction.

I also have a points layer. The points are located at the start of each line (and therefore the end of the preceding line if following in the line direction). The points contain data.

I wish to assign the data associated with each point to the line segment following it, when considering the line direction.

Does anyone know of a way to achieve this utilising QGIS?


texts: Using geometry rather than geography data types in SQL Server Spatial?

texts:
Historically I've mainly worked with simple latitude and longitude coordinates in a SQL Server. I stored them as a geography datatype in SQL Server, and rendered them on a variety of displays (primarily Google Maps).

Recently I've started working with shapefile datasets from various sources, and almost exclusively they are geometry. Further, they use a variety of coordinate reference systems.

It's pretty maddening to work with.

The USGS, for some of their data, use a very rare and hard to ID CRS.
The county of Los Angels says they use the State Plane 5 (without reference a CRS specifically), and of course there are multiple "state plane 5's" to choose from in QGIS.

Are there really some advantages to these particular spatial reference systems that cannot be achieved by sticking with WGS84, considering that the VAST majority of people viewing and consuming will want it in that format?


texts: Why is the performance of ogr methods different on different OS?

texts:
This question actually comes forth from an earlier question of mine
.

A short introduction. I was writing a program that performs near analysis and creates buffers for a set of points. On my Linux Ubuntu machine this turned out to be really slow. Hence the earlier mentioned question of mine.

When running the same code on my OSX it was tremendously fast. At least n comparison to the results I was getting on Ubuntu.

Here are my results

Method     | OSX            | Linux Ubuntu

ogr buffer | 0:00:01.434389 | 0:01:08.384309

sub string | 0:00:19.714432 | 0:00:10.048649

sorted set | 0:00:01.239999 | 0:00:00.600773

So my question really is: why are there such great performance difference with OGR?

Specs Mac OSX
Processor 4x2.5 GHz
Memory 8 GB 1600 MHz

Specs Dell Linux Ubuntu
Processor 8x3.4GHz
Memory 7.8 GB


texts: How to send bbox dynamically via wms request in leaflet?

texts:
I've created a leaflet map and use L.tileLayer.wms to load my wms layer. I dont want to load all data on my map instead  i need to load data based on bbox.I use this line for sending wms request to geoserver :

L.tileLayer
  .wms(elm.GIS_LayersAddressMapServer.trim(), {
    layers: elm.GIS_LayersName,
    format: "image/png",
    transparent: true,
    tiled: true,
    maxZoom: 25,
    bbox: map.getBounds().toBBoxString(),
  })
  .addTo(map);

but in network console i have this :

http://localhost:8080/geoserver/cite/wms
  ?service=WMS
  &request=GetMap
  &layers=cite%3APTL_Tower
  &styles=
  &format=image%2Fpng
  &transparent=true
  &version=1.1.1
  &tiled=true
  &bbox=25.0927734375%2C22.836945920943855%2C78.22265625000001%2C41.50857729743935
  &width=256
  &height=256
  &srs=EPSG%3A3857
  &bbox=5165920.119625352,4226661.916057106,5322463.153553394,4383204.949985149

as you see bbox is sent twice and doesn't show me anything. How can i handle this problem ?


texts: How to find matching pixel values in seperate rasters

texts:
I am looking to scan two separate rasters, capture the pixel information,  then draw a polygon of the matching pixels.

any thoughts? I am at a loss.

Edit:

I have two aerials, both are in different geographic areas (same projection and within about 1 km of eachother). I am looking to determine, if between the images, there are matching pixel grids.


texts: pg_routing source and target column issue

texts:
I have created pgr_createTopology() function on my road data table which accurately assign source and target to destination road table, but I have an issue when I split a road segment in qgis it provide same values of source and targets to both road segments and after that running pgr_createTopology() function is not updating or assign new value for a new split segment which is actually split from old one.

one way which I have found is to null both columns and run pgr_createTopology() again which solves my problem but if I have huge data updating both fields with null again and again is not a good option.

Please provide me solution by which both columns are updating as after I split a road segment in qgis.


texts: Georeferencing Old German Map in ArcMap that uses Ferro Meridian and very uncommon projection

texts:
I am trying to georeference an old map of the German Reich from 1930's. The map is based on the topographische Übersichtskarte des Deutschen Reiches 1:200000.

The coordinates depicted on the map, that I want to use to georeference it, use the Ferro Meridian. I know that the original map uses a conic projection (German: Schnittkegelprojektion) with two parallels at 50° and 53°. Axis of the cone is the same as the axis of the earth. I don't know anything about the central meridian other than that the angle between the central meridian and the corner/border meridians (German: Randmeridiane, do not exactly know the translation) is 0°23'28''5.

The issue is to find a suitable projection, which fits ok and at the same time allows me to input the coordinates in Ferro using ArcMap as the geoereferencing tool.


texts: How do I create attributed line segments between two known points?

texts: I have stream survey data that is spatially referenced every 10th habitat unit by a point.  There are 10 habitat units between each point along the survey and I need to figure out a way to spatially reference each individual habitat unit.  The habitat unit lengths are known from the survey but are not individually referenced.  Does anyone know how to complete this task?


texts: Why are Google Mercator and UTM rotated towards each other?

texts:
Being a self-taught GIS user, I struggled a while before getting the grips on map projections. But now this: because I like to use OpenStreetMap as a background layer for orientation, I have saved my shapefiles all in EPSG:3857 ("Google Mercator").

However, if I want to create a grid of 10×10m² grid cells, I had to convert one of the layers to UTM (zone N32, EPSG:32632) for the spacing to be correct. What surprised me when reactivating the OTF conversion, that the grid had rotated:

I would have expected minor distortion or even scaling, but why is North not identical (straight up) for both projections?

For reference: the excerpt shown is from 48.18714° N, 12.14092° E.


texts:
I am creating a buffer around a PDA point.  Something I have done previously with no problems.  However, this time it is coming up with an error message saying, it has failed to open the Shapefile 'C:\Program Files\QGIS 2,18\bin\temp.shp'.

I've looked into that location and there is no temp shapefile.  Any idea where I can get/re-create one?


texts: Attempted unsupported operation with Raster Catalog in ArcMap?

texts:
I recently installed ArcGIS 10.4 on my Windows 10 machine. I upgraded from ArcGIS 10.3.1.

After upgrading, I noticed when I add a raster catalog to a map and go to the properties, it displays an error, and then the map crashes.

If anyone else has encountered this, any help will do.

I have submitted an error report to ESRI so I will wait to see as well if they respond.


texts: Google maps Javascript API usage statistics does not show

texts:
I am currently using the google Geocoder service from the Google Maps JavaScript API in a few applications.

For months I was receiving usage statistics in my API manager, but for the last 2 months it has gone dark usage screenshot.

I know that the service is still working, and I know that we are receiving around 500 -1000 geocode requests per day across 3 apps... But as we add geocoding functionality to more maps, or if we hit a peak day (we have one or two per year) I am worried that I will come close to my 25,000/day limit.

Anyone know why I lost my usage stats? I have tried the "send feedback" form... This error did occur around the time that I moved my API key from my personal gmail account to our organizations account.


texts: I don't seem to be able to save Layer Styles in QGIS

texts:
I have created and saved layer styles, not recently I admit but have done in the past. For some reason or other I just don't seem to be able to save any new ones in QGIS. No file is created. It does not even save changes to existing ones, even though it says I am replacing it. I surely must be doing something wrong, but don't know what. I can make the changes in the style and that's fine, it just won't save.

I'm running on MS Vista.


texts: Downloading OpenStreetMap and Imagery to ArcGIS Collector for Offline Use?

texts:
I am traveling to Costa Rica and would like to use Collector with my cell phone to record some information in a feature class.  I would like to find some 1m or better imagery to put in a web map that could be synced for offline use.

Does anyone have suggestions for Costa Rica imagery that is publicly available?

Also, is it possible to publish OpenStreetMap as a vector tile package for the same offline use?


texts: Clip (Data Management) tool with basemap

texts:
Can the Clip (Data Management) tool can be used with a basemap?

I tried without success.

If yes, how can it be used?

I am using ArcMap 10.5.1.

I want just cut the image of my property. The sources of images of my property are the basemap "imagery" and aerophotos from my arcgis csw client. With aerophotos (.grid) from other local i could use the Clip tool to cut as i want in my property, but with basemap and client csw i can't.


texts: A number of polygons in .KMZ not visible in Google Earth Pro

texts:
I have downloaded the dataset for the Areas of Outstanding Natural Beauty (England) from various places but some of the areas (at least 5 including Shropshire Hills, Devon, Cotswolds) do not appear in Google Earth Pro. The polygons are listed and GE will zoom to the relevant polygon but is not visible.

Data can be found at:
https://data.gov.uk/dataset/areas-of-outstanding-natural-beauty-england3/resource/23d670a5-9c42-4cca-a2db-1c1d910d4b8f

Incidently when I proceed to fill the layer it also fails to work on other polygons - (e.g. Quantock Hills, Chilterns)


texts:
I try to get the QGIS Server running and just try to establish a first WMS.
I tried to follow
http://hub.qgis.org/projects/quantum-gis/wiki/QGIS_Server_Tutorial#Windows
and
http://anitagraser.com/2012/04/06/qgis-server-on-windows7-step-by-step/

I already have problems with the installation. All manuals seem to be outdated. I downloaded the osgeo4w file and it starts with installing QGIS 2.5 I would have thought that it installs QGIS 2.4. But what the heck. Then all manuals refer to the folder
c:/osgeo4w/apache/conf/httpd.conf
But the installation via osgeo4w does not include an installation of apache anymore.

So no I do not know how to go on. Do I have to install apache in my own or is there another solution without the apache?

I just tried to install apache to c:/osgeo4w/apache/ with the port 81. I get the message "It works!" in my Browser.
I wanted to copy the files

libeay32.dll and ssleay32.dll
from C:\OSGeo4W\apache\bin
to C:\OSGeo4W\apps\qgis\bin.

As mentioned in the first "How to" but these files does not exist in this directory and they do not exist in any directory of apache. These two files are in:
C:\OSGeo4W64\bin

I tries to copy the in both first mentioned directories and tried to start the adjusted capabilities request mentioned in the "How to", but the the file is not found on the server

localhost:81/qgis/bin/qgis_mapserv.fcgi.exe?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities

What can I do to make an WMS running via QGIS?


texts: Pixel Coordinate to latitude and longitude conversion

texts: I have a map of part of Brazil. I want to find the latitude and longitude of some marked points on the map. Is there any way to find latitude and longitude of a point from its pixel coordinates given the map is using Mercator projection? I found a similar question in Georeferencing old map in QGIS? but there lat-lon grid was provided but here in my question nothing is given .


texts: How to set focus on ArcMap map window using arcobjects and c#.net

texts:
I have developed a custom dockable window with textboxes, comboboxes and a button for use in ArcMap. After the user clicks on the button in the dockable window a function is invoked. After the function has finished, the focus is still on the dockable window. That’s why the user has to click into the map window, to set the focus to the ArcMap map window.

Now the question: Is it possible to set the focus to the ArcMap map window programmatically? Thus, the user could immediately use keyboard shortcuts to navigate or edit the map and does not have to click into the map before. I would implement such a functionality into the function invoked by the button in the dockable window.

Until now I haven’t found any arcobjects interface or property which does the job. And I have no idea how I could get the handle of the ArcMap map window.

Thank you


texts: Lagrange Multiplier in lm.LMtests (spdep package for R)

texts:
Are the Lagrange values comparable in lm.LMtests?

If both RLMerr and RLMlag are significant, would the test with the largest Lagrange value be considered 'better' at explaining the spatial relationship for the given model? Or 'better' at explaining which type of spatial autocorrelation is most prevalent in the model?

This is my first question, so apologies for probably not following the necessary protocol/formatting/salutation OR if this belongs in some other silo of StackExchange.


texts:
I am trying to implement Dijkstra's algorithm, my query is working fine if I set direction parameter to false. I have to consider the direction of travel while calculating route from point A to point B. This algorithm determines the direction by the order of the coordinates in the linestring, this is where I am facing problem. If I am traveling from north to south and the linestring coordinates of that route are digitized from south to north then the algorithm is not returning the correct route. I am using geometry length as cost and reverse_cost.

So, my question is how can I solve this problem? I have another column in my table containing Boolean value if the road is one way or not, can I use this column somehow to find the correct route?


texts: Change suffix label on decimal grid in print composer

texts: I want to change my suffix label on decimal grid/measured grid coordinates for a map canvas in the print composer from "N" or "E" to "m.N" and "m.E". How can I achieve this?


texts:
I'm very new to python so please bear with me. I am trying to run a model that exports a .wm feature to a geodatabase feature class. The .wm file extension is for software that is used to update electrical utility data. The model has to be run in an .mxd because the .mxd contains an add-on that allows arcmap to read the .wm files. I'm having a problem getting the script to run successfully because it is not recognizing the .wm file even though the .mxd is open with the add-on enabled. Is there a way to prompt the model to run from the mxd rather than separately? Below is my script. Any insights?

Script:

import os
import arcpy
mxd = (r"D:\Milsoft\WindMilMap\WindmilMapExport\WindmilMapExport.mxd")
os.startfile(mxd)

import time
time.sleep(50)

import arcpy
arcpy.ImportToolbox(r"D:\Milsoft\WindMilMap\WindmilMapExport\WindmilAutomation.tbx", "WindmilAutomation")
arcpy.ElectricWMtoGDB_WindmilAutomation()

Traceback:

Traceback (most recent call last):
  File "D:\Milsoft\WindMilMap\WindMilMapExport\WindmillExportAutomation.py", line 11, in <module>
    arcpy.ElectricWMtoGDB_WindMilAutomation()
  File "D:\Milsoft\WindMilMap\WindMilMapExport\WindmilAutomation.tbx", line 21, in ElectricWMtoGDB
       F e a t u r e   L a y e r
ExecuteError: Failed to execute. Parameters not valid.
The value does not exist.
WARNING 000258: Output D:\Milsoft\WindMilMap\WindMilMapExport\Electric.gdb already exists
ERROR 000732: Input Features: Dataset 'D:\Milsoft\WindMilMap\WindMilMapExport\WindMilMap LLM E-CMS.wm\Capacitor' does not exist or is not supported
Failed to execute (ElectricWMtoGDB).


texts:
I'm trying to find islands (type not water, completely surrounded by water) in a polygon table. The polygons are not overlapping.

Data:

CREATE TABLE islands (
   ID serial PRIMARY KEY,
   type VARCHAR (10),
   geom geometry);
INSERT INTO islands(type, geom)
   VALUES
   ('water',st_geomfromewkt('SRID=3006;POLYGON((510796 6429408,513939 6420063,526040 6421096,527763 6430183,517341 6436643,510796 6429408),(516954 6429623,516652 6425747,518203 6426436,519667 6428073,516954 6429623))')),
   ('land',st_geomfromewkt('SRID=3006;POLYGON((516954 6429623,516652 6425747,518203 6426436,519667 6428073,516954 6429623))')),
   ('land',st_geomfromewkt('SRID=3006;POLYGON((510796 6429408,517341 6436643,512561 6436427,508384 6434619,510796 6429408),(511011 6434231,512475 6435652,513207 6433025,511011 6434231))')),
   ('land',st_geomfromewkt('SRID=3006;POLYGON((511011 6434231,512475 6435652,513207 6433025,511011 6434231))'));

Only id 2 should be selected:

I have tried query below based on Determine if a polygon is not enclosed by other polygons, but it is not selecting anything.

SELECT a.id
FROM islands a
INNER JOIN islands b ON (ST_Intersects(a.geom, b.geom) AND a.id != b.id)
where a.type = 'water' and b.type <> 'water'
GROUP BY a.id
HAVING 1e-6 >
  abs(ST_Length(ST_ExteriorRing(a.geom)) -
  sum(ST_Length(ST_Intersection(ST_Exteriorring(a.geom), st_interiorringn(b.geom,1)))))


texts: Calculating proportion of pixel containing vegetation (from sentinel 2 image stack)

texts:
My task is to perform a supervised classification of landcover across a huge (20,000+ hectare) property. I need to classify into three classes: forest, regenerating trees, and soil. We use s
Sentinel-2 imagery for this task.

However, both the randomforest classifier and maximum likelihood classifier get confused as some regenerating trees are as small as 1m2 and the image is 5m2. so they're surrounded by soil and we can't classify them accurately.

Is there a way to calculate the amount of tree in each pixel?

My current thought process is that we could make training samples that are as accurate as possible (and close as possible to spectrally pure). Then we could set a really narrow threshhold in the classification so only pixels which we are certain are tree/ground are classified. Then with the remaining unclassified pixels  we perform some form of calculation based on colour to find how much vegetation could be in the remaining pixels.


texts: Arcmap 10 Dot Net SDK Installation Issues

texts:
We have had Arcmap 10 installed on our machines for about a month now and all seems to be well.  I ran into an issue yesterday when I tried to install the Developer SDK for .NET.  An error message shows up stating that I need to have a supported version of Visual Studio installed before the install can continue.

We do have Microsoft Visual Studio 2008 Sp1 Professional installed as well as .NET 3.5 framework.  Visual Studio and Arcmap both work just fine.  Has anyone else ran into this problem or know of a solution?  I have put in a help desk request with ESRI, and will post here if I get an answer from them.


texts:
I'm running a i3 3,1 GHz Machine, 4 GB of RAM with an onboard Graphic Chip.

I build a small GIS with Mapbender3, Geoserver based on PostGIS data.
It's running pretty fluid from what i can see. Now i want to insert high-resolution Satellite Images i have on my hard disk instead of Loading a WMS providing those images.

Do you think my hardware will go to its limits in that case ?
I have very small GIS-Data, like 100 datasets.

I told my Computer merchant to make me an offer for a PC that can handle this System and he made an an offer with following components:

GB GA-Z97P-D3 Mainboard
intel Xeon E3 1231V3 3,4 Ghz
16 GB RAM
Sapphire FirePro W5100 4096MB

To me it looks like, this beast would be heavily overpowered ?
is such a Graphic card necessary to handle two dimensional data ?


texts:
I'm curious if, and how, people use EPSG defined water level vertical datums such as Mean Low Water.

There is specific guidance to not use the OGP EPSG definition:

  Remarks: Users are advised to not use this generic vertical datum but to define explicit realisations of MLW by specifying location and epoch, for instance "MLW at xxx during yyyy-yyyy".

Which is a blanket statement across all water level vertical datums sourced from the IHO Dictionary, UK Hydrographic Office, and Canadian Hydrographic Service.

Is there a common convention, especially with WKT, to express a datum as a realization of an EPSG defined datum?


texts: Is there a way to search raster files with a bounding box using GDAL in C#?

texts:
So far I have tried this
        GdalConfiguration.ConfigureGdal();
        var outputFile = "C:/Users/abcd/RPF/output.vrt";
        var rasterFiles = Directory.GetFiles("C:/Users/abcd/Desktop/RPF/xxx01/").Where(x=>x.EndsWith(".ON1")).ToList();
        var vrtOptions = new GDALBuildVRTOptions(new[] { "-overwrite", "te 1.91241,64.9814,3.1241,101.101", @"input_file_list C:\Users\abcd\RPF\list_of_raster.txt" });
        var vrtDataset = Gdal.wrapper_GDALBuildVRT_names(vrtFile, aeroFiles.ToArray(), vrtOptions, null, null);
        vrtDataset.FlushCache();

This creates a .vrt file but with all the files present in the directory, but I know the bounding box co-ordinates matches only some files. How can I do this?


texts: Is there an easy way to implement INSPIRE model in PostGIS?

texts:
I need to implement INSPIRE (Infrastructure for Spatial Information in the European Community) model in PostGIS.

I searched everywhere but couldn't find SQL file, or any other implementation that I could download. Is there any way to quickly auto generate empty tables with all relations, or do I have to manually generate 400+ tables?


texts: Restricting viewshed to less than 360 degrees?

texts:
I'm using QGIS with the visibility analysis plugin.
I want to do a viewshed in only one direction. Is it possible to restrict the analysis rather than doing a 360 degree analysis?


texts: Setting query loop in ModelBuilder/ArcPy?

texts:
I am looking to apply a Select By Location protocol to my data, but need to first query individual IDs in my target and source layer (the same ID in each table, one at a time). How do I loop through a query? I have far too many IDs to manually set a query in each table then perform the selection.

Is there a way to create a looping model (with parameters) or script that will systematically select each (matching) ID in my target and source layers then perform the tools I specify?

I have used the Iterate Feature Selection tool in ModelBuilder to select by ID field in my target layer, but ModelBuilder will not let me have two iterators in the model (even if I could, there is no way to guarantee the iterators are selecting the same IDs at the same time). I have tried linking the Value output of the iterator to a Select by Attribute tool for my source layer but the value can only be set as the whole expression, not part of it.
So I am nearly there.

After extensive discussion with Roy this is the code I've come up with. Unfortunately, while no errors pop up when I run it, nothing actually happens. The process is marked as completely immediately after I run it and nothing has changed. Any thoughts?

rows = arcpy.SearchCursor("hmb_logbook_2009")
for row in rows:
    whereClause = '"trip_id" = ' + str(row.trip_id)
    arcpy.MakeFeatureLayer_management("hmb_logbook_2009", "currentLines", whereClause)
    arcpy.SelectLayerByLocation_management("hmb2009", "WITHIN_A_DISTANCE", "currentLines", "1 Miles", "NEW_SELECTION")
    arcpy.SelectLayerByAttribute_managemet("hmb2009", "SUBSET_SELECTION", ' "trip_id" = ' + str(row.trip_id))
    arcpy.CalculateField_management("hmb2009", "match_activity", 1)
    arcpy.SelectLayerByAttribute_managemet("hmb2009", "CLEAR_SELECTION")
    arcpy.Delete_management("currentLines")
del row, rows


texts:
In Pablo's answer to How to open a GeoJSON file with mixed geometries in QGIS?

He says he 'imported each type at a time' and I do not know how to do that.

I am working with a GeoJSON file that has a mixture of points and polygons, but
points are first, so it only shows points.

How do I get QGIS to show the polygons instead of the points?


texts: QGIS not launching after new installation on SuSE Linux

texts:
I have a brandnew installation of SuSe Linux leap 15.2, and installed the latest version of QGis (3.18).
Installation with Yast runs fine, no errors.
However, QGIS does not launch from the start menu. Launching from the console gives an error message indicating a problem with the qgis library:
qgis: symbol lookup error: /usr/lib64/libqgis_core.so.3.18.2: undefined symbol: GEOSCoordSeq_setXY_r
I tried also replacing the installation with the QGIS ltr version (3.16) but that gives the same type of error
Any tips how to resolve it?


texts:
I have evaluated the Getis-Ord Gi*-spatial statistic for each point in my map. Like this, each point is assigned a Gi-score and due to this a certain color, as you can see in the image below.

What I would like to do now is to group points of the same color and make them appear as an area of that color, similar to the second image.


texts: QGIS rule-based categorization with multiple variables

texts:
There are a few questions which sort of touch on this, but I haven't found an answer that is easy yet.  I have a database with several columns of metadata.  The first column relates to a unique ID I have assigned to point which can be shared among many points.  The second column is the original data source of the point which again can be shared by multiple points.

I would like to plot this data so that the unique ID is categorized by colour and data source is categorized by shape.  This can be accomplished using rule-based categorization, however I have 70,000 data points with 50 sources and 1500 unique ID's; setting up these rules by hand is not feasible for this project.

Is there a simpler way to do this?


texts: Tiger geocoder is not finding results from manually loaded tables

texts:
I'm trying to test the TIGER geocoder with the address provided in the documentation.

This is the sample address I'm testing:75 State Street, Boston MA 02109

SELECT g.rating, ST_X(g.geomout) As lon, ST_Y(g.geomout) As lat,
    (addy).address As stno, (addy).streetname As street,
    (addy).streettypeabbrev As styp, (addy).location As city, (addy).stateabbrev As st,(addy).zip
    FROM geocode('75 State Street, Boston MA 02109') As g;

But I'm not getting any match for this query:

I have installed PostgreeSQL, PostGIS and the requerired extensions:

CREATE EXTENSION postgis;
CREATE EXTENSION fuzzystrmatch;
--this one is optional if you want to use the rules based standardizer (pagc_normalize_address)
CREATE EXTENSION address_standardizer;
CREATE EXTENSION postgis_tiger_geocoder;

Also I have also downloaded the data corresponding to that address from the census site, which is correspondes to the state of MA (suffolk_county)

I have loaded the data into a table using the following comand:

ogr2ogr -f PostgreSQL PG:dbname=test1 -progress -nlt PROMOTE_TO_MULTI /GIS/data/ma_suffolk_county/tl_2016_25025_edges.shp -nln ma_suffolk_county

The table was succesfully created:

\d+ ma_suffolk_county


texts:
I would need to create a complex legend divided into many sub-sections that does not correspond to different layers but to different rule-based styling objects. It seems to me that this possibility is not supported by the legend items in print composer, how can I overcome this problem?

p.s. the images refer to a part (some layers) of the entire legend


texts: How to merge large number of vector features with QGIS SingleToMultipart?

texts:
My goal: Merging all line segments in a vector to a single feature.

My approach: QGIS module singlepartstomultipart

My problem: For larger numbers of features within the line vector, the module has a processing time of several hours or even makes QGIS crash.

My code:

  processing.runalg('qgis:singlepartstomultipart', {"INPUT":
  vect_in, "FIELD": "DN", "OUTPUT": vect_out})

My question: Is there a different way for solving this?

My suggestion: I tried to run singlepartstomultipart on sets of 100 or 1000 features to then merge these again. However, this seems to be quite a detour and running singlepartstomultipart didn't work on the merged sets (see #? in the code below. This code is the implementation of my suggested detour solution):

# Open line layer
Lines = QgsVectorLayer(Lines, Lines, "ogr")

# Define layer extent
extent = Lines.extent()
xmin = extent.xMinimum()
xmax = extent.xMaximum()
ymin = extent.yMinimum()
ymax = extent.yMaximum()

# Initialize counter
j = 0

# Merge line segments in sets of 1000 features
i = int(round(float(Lines.featureCount()) / 1000))

while j < i:
    # Create temporary layer
    vl = QgsVectorLayer("LineString?crs=EPSG:4326", "temp_layer", "memory")
    vl.setCrs(Lines.crs())
    pr = vl.dataProvider()
    vl.startEditing()

    # Select block of 1000 features to be merged
    start = j * 1 + j * 1000
    end = (j + 1) * 1000

    # Copy features to temporary layer
    iter = Lines.getFeatures()
    for feature in iter:
        if feature.id() > start and feature.id() < end:
            pr.addFeatures([feature])
            vl.commitChanges()

    # Add layer to registry
    # Note: Layer needs to be added to the canvas in order to be usable in processing.runalg below (#unclear)
    QgsMapLayerRegistry.instance().addMapLayer(vl)

    # Merge line segments to one feature
    multipart = os.getcwd() + r"\1a_multipart" + str(j) +".shp"
    if not os.path.isfile(multipart):
        processing.runalg('qgis:singlepartstomultipart',
                          {"INPUT": vl, "FIELD": "DN", "OUTPUT": multipart})

    # Remove layer from registry
    QgsMapLayerRegistry.instance().removeMapLayer(vl.id())
    j += 1

# Merge created files
list = []
for file in os.listdir(data_dir):
    if '1a_multipart' in file and '.shp' in file:
        temp_lines = QgsVectorLayer(file, 'temp_lines', "ogr")
        list.append(file)
multipart_merge = os.getcwd() + r"\1b_multipart.shp"
if not os.path.isfile(multipart_merge):
    processing.runalg('qgis:mergevectorlayers', list, multipart_merge)

# Merge features in merged shapefiles
final_merge = QgsVectorLayer(multipart_merge,multipart_merge, "ogr")
multipart = os.getcwd() + r"\1c_multipart.shp"
if not os.path.isfile(multipart):

    #? This doesn't work so far
    processing.runalg('qgis:singlepartstomultipart',
                      {"INPUT": final_merge, "FIELD": "FID", "OUTPUT": multipart})


texts: How to remove black background off raster images?

texts:
I am having trouble with QGIS 3.4. I have just added a raster layer of a JPEG image of an island I have stitched together, however on QGIS it is displayed with a black background around the edge. I would like to remove this, I understand I might have to reclass the data?

Is it showing up as black as there is no data?