texts: Understanding Join Data vs Spatial Join?

texts:
I am working on a semester project for my GIS Programming and Customization course. I built a model that receives updated crime point data from the City of Denver as an input.

The sequence of tools that I have strung together in the model are as follows: Select Layer by Attributes (because I only want to use crime points that represent reports of domestic violence), Create Feature Class, Add Field (this is a numeric field named "Count"), and Calculate Field (populates all records with the value of 1).

Now I have a new feature class that consists of points where domestic violence crimes have been reported and each record is populated with the value of 1 in the "Count" field. I want to join this data to a separate polygon feature class (Census Neighborhood Demographics), so that my resulting polygon feature class ends up with a "Sum Count" field that represents the total number of domestic violence crimes that have been reported in that neighborhood.

I tried including the Spatial Join tool in my model in order to complete this step, but when I run the model (or try using the tool manually), my output polygon feature class still shows every record populated with the value of 1, rather than the sum of crime points.

I have done this whole process manually before, but by right-clicking on my target layer and choosing to join it to another layer based on spatial location. I noticed that this operation is called "Join Data". The dialog window is very different than that of the Spatial Join tool. I have tried searching for a tool called "Join Data", with no success.

Does anyone have any suggestion for how I can make my model do what I was able to do manually by right-clicking on my target layer and opening the "Join Data" tool?

I am very open to using a snippet of Python script, but I am new to the programming language and am not sure how to accomplish my goal with code.


texts: Calculate density of point feature class using field calculator? (not pop)

texts:
I'm just wondering what would be the best method to calculate density in the field calculator for a point shapefile (ArcMap). This wouldn't be based on a population value. I would need to derive a value for each point based on the proximity of surrounding points that I can plug into each point feature. Obviously the more points surrounding the subject point the higher the value. How would I do this?
I have the spatial analyst tool
I don't want to create a raster. I want to be able to symbolize each point based on their density value (stretched)


texts:
I've created a text field, width 255, called BUFFER in the feature class "ForestRoads_Buffer".  I want to calculate the BUFFER values of all the features to "Y".  The code below doesn't do this and I don't understand what the expression is supposed to look like.

arcpy.CalculateField_management("ForestRoads_Buffer", "BUFFER", "Y", "VB")


texts: Esri Python Toolbox(pyt) corruption in ArcToolbox.dat

texts:
Symptoms:

Open ArcCatalog and all tools within Python Toolbox(pyt) have a red-x error next to them claiming that "Arcpy module has no attribute Parameter"
Close ArcCatalog
Delete ArcToolbox.dat file in ~\AppData\Roaming\ESRI\Desktop10.2\ArcToolbox
Open ArcCatalog, everything is fine with tools
Run tools successfully
Close ArcCatalog
GOTO Symptom 1 .....

Seems like something in the geoprocessing results history is going wrong.

Any ideas for fixes?

I also want to maintain geoprocessing result history.


texts:
I am writing a JSON file, as string input, using python, into a .json file in the file system.
However I am getting the following error sometimes and I don't know why.
Any ideas?

  UnicodeEncodeError: 'ascii' codec can't encode character u'\xe4' in
      position 1920: ordinal not in range(128)

The json file created is used as input in the "JSON To Features" geoprocessing tools and it fails.

# Script arguments
inputFeatures_json = arcpy.GetParameterAsText(0)
#write to file
jsonFileName = 'WebMap_{}.json'.format(str(uuid.uuid1()))
JSON_File = os.path.join(arcpy.env.scratchFolder, jsonFileName)
f = open(JSON_File , 'w')
f.write(inputFeatures_json)


texts: gdaldem slope....odd values, almost all 89.98+/- (? - why?)

texts:
I'm passing an SRTM1 3601x3601 file to "gdaldem slope".....I've queried the values in Python, they are reasonable (min 222, max 424) for an area of SW Iowa.  Took a bit to get the syntax right, but file is being created all right.  Oddly though, almost all cells have a value from say 89.95 to 89.99, with a few having zeros (and the borders, no surprise / no problem, 9999).  Looking at the 0's, it is when the elevation is flat, say for a water body.

Default slope is supposed to be degrees; anyone have any idea why I am getting the ~89.9x values?


texts: ArcPy TableToNumPyArray where_clause issue

texts:
I am trying to run this code in the Python terminal in ArcGIS, but everytime it throws an error and I can't figure it out. I think it probably has something do with how I construct the SQL where clause because when I remove that it works fine.

import arcpy
import numpy as np

##VARIABLES TO SET
outputFile = "F:\\output\\data.csv"
tableName = "myTable"
pointsName = "folder\\myPoints"

pointsArr = arcpy.da.FeatureClassToNumPyArray(pointsName, ('id', 'KEY'))
keyStr = ", ".join(map(str, pointsArr['KEY']))
sql = "{0} IN ({1})".format(arcpy.AddFieldDelimiters(tableName, 'key'), keyStr)
tableArr = arcpy.da.TableToNumPyArray(tableName, ('key', 'column'), where_clause=sql)
numpy.savetxt(outputFile, tableArr, delimiter=",", fmt='%s')

This is the error I get when I try to run it in ArcGIS.

Traceback (most recent call last):
  File "<string>", line 12, in <module>
RuntimeError: The database was not found.

Can any help me figure out how to get this working? Basically what I'm attempting to do is get the keys from my 'points' layer and then use those keys to query a standalone table.


texts:
We have a quite large ArcSDE (~120 versioned feature classes & tables, 1M-3M records) which has Oracle or SQL underlying DBMS. This SDE is exposed through a MapService and we use ArcGIS JS API to display layers.

Now we'd like to have full-text search against set of feature classes and tables in the SDE. ESRI manuals say we should set up XML columns to enable search through documents. But I really don't understand how it helps.

So the questions are:

Is ArcSDE capable of full-text searching?
And if yes, is there any services in ArcGIS Server or whatever which can help us with this?

Firstly, I'd like to say that I little bit messed up terms. I thought that ArcSDE==SDE==geodatabase. However It's different things and I actually meant geodatabase.

Secondly, here's the solution we've came up with: There's  Find task available on ArcGIS Map Service out of box, which does text search against given layers. It performs unexpectetly fast and we're going to use it. If performance is unsatisfactory we'll enable Lucene to cache queries. Hope this solution will work. I'll report the result here as soon as we implement this approach.


texts: Assign points to another point layer's value

texts:
ArcMap 10.2 with extensions ArcScan and Spatial Analyst.
I have two point layers (A and B). A has about 40 points and each point has a value (max. number of people it can hold). B has about 1000 points (every point is a person).
I want to assign B (people) to A's max. number so that I can tell if there is enough capacity in all A-points for all the people on my map.
How do I do this, do I need to convert it to raster data or can I just do an analysis with the point data?
(I haven't been using GIS for a few years so I'm a bit rusty. I hope this question isn't too stupid.)


texts: How to handle two onclick events that not both are fired?

texts:
I make a web application with ArcGIS API for Javascript (version 3.15) to access a geoprocessing service. The service is executed on an onclick event by the user, because the clicked position is used as an input parameter:

map.on("click", startGPService);

Then the result is displayed and I have an info window enabled, which also opens when clicked, although not specified directly by onclick:

var template = new InfoTemplate();
template.setTitle(getTitleInformation);
resultLayer.setInfoTemplate(template);

The problem is that whenever the user clicks on an info window, also the onclick event of the map is fired and the geoprocessing starts again. It would be great, if the onclick event could recognize, if the info window is clicked or somewhere on the map and only then fire.

Any solution how to deal with that?

This question is linked to this topic: ArcGIS API - Creating Popups for Geoprocessing Result


texts: Can we do replication and versioning in User schema geodatabases version with sde default version

texts: As both are administrators to their own schemas, is there any possibility to post and reconcile data across schemas?


texts: Looping through table to get single intersection from N>2 geometries using PostGIS?

texts:
It's really clear how to apply ST_Intersection when you have only two geometries. But I am struggling with the situation where you have N>2 shapes and you want to compute a single one that represents intersection of all of them as a single geometry.

With the help of ST_Intersection: Intersection of all geometries in a table I can get intersections for all pairs of geometries in a single table. But it is not exactly what I am trying to achieve though.

Let me try to illustrate the problem. What I have is this (simplified to only three geometries):

But what I am after is this:

I have over 5000 shapes in the table, so obviously I can't specify them explicitly.

I imagine a loop could be a solution, something like:

result = union(geometries)
for geom in geometries:
    result = intersection(result, geom)

How would you achieve this with PostgreSQL/PostGIS?


texts: Numpy TypeError when turning ndarray into structured array

texts:
I want to turn my ndarray(670911L, 39L) into a structured array using a view:

c = arr.view(dtype=[(n, ‘float64’) for n in names]).reshape(len(arr))

where names is a list of 39 strings. But I am getting a TypeError: data type not understood.

What am I doing wrong?


texts: Adding existing style while publishing a new layer in GeoServer using REST API

texts:
I have a function which publishes a layer in GeoServer with the REST API.
Now I want to add an existing style in the layer while publishing it.

Anyone has any idea where I write the style name or call the 'style.sld' while publish the layer?

const publish = function (layerName) {
return new Promise((resolve, reject) => {
    axios({
        method: 'POST',
        url: `http://${GEOSERVER_IP}:${GEOSERVER_PORT}/geoserver/rest/workspaces/workspaceName/datastores/datastoreName/featuretypes.xml`,
        auth: {
            username: `${GEOSERVER_ADMIN_USER}`,
            password: `${GEOSERVER_ADMIN_PASS}`
        },
        headers: {
            'Content-Type': 'text/xml'
        },
        data: `<featureType>
                <name>${layerName}</name>
                <srs>EPSG:4326</srs>
            </featureType>`
    })
        .then(response => {
            resolve({
                ok: 1,
                message: `Layer has been published.`
            })
        })
        .catch(err => {
            reject({
                ok: 0,
                message: err.response.data,
            })
        })

})
}


texts: Using update cursor on a field with a domain

texts:
I have been replacing field calculations with python update cursors in many of my work flows and today I ran into a problem.  The fields being updated have an attribute domain.  so now, when I "calculate" the field based on the domain value, the field is not populated with the desired description.  Instead it is populated with the value used in the calculation.  Please see the portion of code for the update cursor.

# Get linear units
(path, name) = os.path.split(featureClass)
spatial_ref = arcpy.Describe(path).spatialReference
units = spatial_ref.linearUnitName
if units == "Meter":
    lengthUOM = "Metre"

# Calculate feature class fields
rows = arcpy.UpdateCursor(featureClass)
for row in rows:
    length = row.getValue("Shape_Length")
    row.setValue("featureLength", length)
    row.setValue("UnitOfMeasure", lengthUOM) #this field has a domain
    rows.updateRow(row)
del row, rows

When I use field calculator (in an edit session) and set the value to "Metre" the field will populate with the description from the domain (what I would like).  The above code, however, sets values in the field to "Metre" instead of the description.  I have read this question, but it does not address update cursors.

Has anyone had any success with update cursors on fields with domains?  If so, could you give me some insight into where lies my mistake?  Thanks in advance for any advice.


texts: Cropping raster and merge it to overlapping DEM?

texts:
I am working with an existing DEM of a river and its surrounding area in QGIS (actually a virtual raster from four separate DEM TIF tiles).

My goal is to add a bridge's profile to the DEM. I have a point cloud (.las) of the area, and from that, have extracted the bridge's area using Cloud Compare using the clipping tool:

When I rasterize the bridge (in Cloud Compare), it also creates an elevation model around the bridge, and appears to assign it at sea-level. When I bring the bridge raster and the DEM into QGIS the bridge raster includes the sea-level generated elevation (outlined in red):

So my question is, how do I extract just the bridge and insert it into the existing DEM?

I mean, I guess I could just rebuild the entire DEM since I have the LAS files and include the bridge elevation, but this doesn't seem like it should be that hard. (having said that, I have spent the morning banging my head against this).


texts: Can I create a Styled Map for only one section of Google Maps?

texts: I have been reading everything I can to find out how this can be done.  It looks like the styled map applies only to the whole map, and I can't find a way to only apply it to a certain polygon region.


texts: Populating dictionary using ArcPy creates only 1 entry?

texts:
I am new to arcpy, and I am trying the simple operation of creating a dictionary that will pair my fcs in a gdb with their corresponding coordinate system.
The operation itself works, but it only returns ONE (random) entry in my dictionary, instead of all fcs in the gdb, and I can't figure out why.

My code reads as follows:

# Import arcpy module
import arcpy
# define source of data
arcpy.env.workspace = "R:/02_Vektor/ALKIS.gdb"
arcpy.env.overwriteOutput = True
# list fcs
fcs = arcpy.ListFeatureClasses()
# create a dictionary that couples shapes with their respective COS
for fc in fcs:
    COS = arcpy.Describe(fc).spatialReference.name
    MK_dict = {
        fc : COS
        }
# check output
print MK_dict


texts: PyQGIS - can't changeAttributeValue - returns False

texts:
Part of my QGIS plugin is responsible for appending a new column to a vector layer and fill this column with values calculated earlier.

I'm able to append the column using

layerPoly = QgsMapLayerRegistry.instance().mapLayersByName(PolyName)[0]
layerPoly.dataProvider().addAttributes([QgsField(newName, QVariant.Double, "double", 10, 2)])
layerPoly.updateFields()

So far it works, creating a new field of a double type, length 10 and precision 2.

Now starts the part of updating attributes of each feature:

columnid = layerPoly.fieldNameIndex(PolyName)

layerPoly.startEditing()

for feature in range(len(self.polygons)):
    # self.polygons is a list of features from layerPoly

    feat = self.polygons[feature]
    fid = feat.id()

    # self.y is a list of values that has length same as self.polygons
    item = round(self.y[feature],2)

    layerPoly.changeAttributeValue(fid, columnid, item)

layerPoly.commitChanges()

print type(fid) returns type: long instead of QgsFeatureId as I think it should?
print type(columnid) returns type: int
print type(item) returns type: float

What causes the changeAttributeValue to be False? I thought it's fid's type not QgsFeatureID but int, but compared to other codes I have seen on the Internet, this fid, feat.id() part is quite the same, so it shouldn't be an issue.

I'm using QGIS 2.8 and PyQt4==4.11.4.


texts: Error message 1 Failed to repack shapefile in QGIS 3.10.3 when saving changes to layer

texts:
I am self taught using QGIS 3.10.3 on Windows 10.
Over time I have a project dealing with 10 Parishes with maps to show land use and ownership using 20 layers that take data from individual CSV files, one for each parish.
I have recently redone the CSV files and relinked the files to the relevant layer. When checking the layers I find some ploygons have vanished. The data is still in the Atribute Table but the polygon does not show.
If I redraw the polygon it still does not show but data is added to the Atribute file. When saving the layer I am now getting a message saying

Error 1 "Failed to repack the Layer".

I have modified the polygon and saved it and I dont get the message but if I try again I do. I keep saving the project file but I am having to keep adding polygons to get he map to draw and then delete the extra data from the atrbute table.
I have reloaded QGIS.
I have also  tried removing a laver, say land use,  and worked on the other owners layer for the Parish When completed I duplicate the Owners layer and rename itan land use but the same issue arises with "Failed to rrepack" message. Advice please


texts: Transforming from NGVD29 to EGM96/WGS84?

texts:
I have this USGS reservoir elevation in NGVD29 datum and I want to transform it to WGS84 ellipsoid (SRTM 1 ArcSec DEM).

How can I transform this point data (NGVD29) to EGM 96 or WGS84 ellipsoid so that later I can calculate the associated volume of water?


texts: How to change timeslider step in CartoDB/torque?

texts: I am using CartoDB and I have a torque map that I am using for a course.  The visualization only covers three months of time.  I would like the timeslider to show each day, but right now it shows every two days.  I have tried manipulating the duration and steps, but without success.  As a non-programmer, is there any why for me get the visualization and timeslider to display every day rather than every two days?


texts: Why does my RAM cap out around 150 MB when doing the Distance Matrix analysis in QGIS?

texts:
I am currently trying to determine the 2 nearest points for a database of 8 million people against 425 locations.

I have 2 csv files uploaded into QGIS one with a unique number for 8 million people with a latitude and longitude. The second csv has 425 locations with a latitude and longitude. I then start the distance matrix analysis under the QGIS Geoalgorithms section of the processing tools feature.

Once I set it up and tell the analysis to run it begins. However, it takes about 4 hours to process which is not bad when considering the amount of calculations I am asking, but when I open the task manager QGIS is only using 20% of CPU and 150 MB of RAM.

To test it further I opened up to instances of QGIS running the same distance matrix and both are running at about 150 MB.

Is there a known limit of RAM use for vector analysis where it is reading from a CSV and writing to a new csv output?

I am currently using QGIS 2.14.15-Essen on Windows 8 with an i7 and 32 GB of RAM.


texts: Adding count of features in categories to legend using ArcGIS for Desktop?

texts: I am using ArcGIS 10.1 for Desktop and I would like to know How to add the Count (which is visible on the Symbology tab of the Layer Properties dialog) to the legend?


texts: Creating Lines in Bounding Box at specified distance apart?

texts:
How would I automate creating lines a specified distance apart with a square shaped feature class, oriented East to West?

I'm using QGIS.

I tried searching, however I haven't been able to find something specific to what I'm looking for.


texts: Populating new field in shapefile with values from list?

texts:
I'm looking to populate a newly created field in a shapefile that from a list using basic python and osgeo tools. The values are integers and in a python list that is the length of the number of features in the shapefile. The goal is to have a set of unique numbers for grabbing random samples of features. Here is example.

from osgeo import ogr, osr
import os, sys
import random

driver = ogr.GetDriverByName ("ESRI Shapefile")
sampleshapefile = ogr.Open('shapefile.shp', 1)
sample = sampleshapefile.GetLayer()
new_field = ogr.FieldDefn('randnum', ogr.OFTInteger)
f40.CreateField(new_field)
ran = random.sample(range(1, len(sample)+1), len(sample))

I've seen the examples such as Programmatically populating shapefile field with polygon areas in OGR?. Yet in that example the for loop on each feature fills the field with the result of a geometry calculation on the feature. I'm looking to have the first value in the list "ran" be value in the "randnum" field of the first feature, second list value for the second feature and so on.


texts: OrfeoToolbox SVM classifier module doesn't produce a model for final Prediction

texts:
I'm having a struggle with SVM classification in QGIS OrfeoToolbox.

I used to work with RandomForest in R, but now I'm trying to figure out how the SVM algorithm works. I guess there would also be similarity between those packages(another topic).

Because I've no scripts nor knowledge about "e1071-SVM-Kernel" package, I've found a plugin within OrfeoToolbox; it simply requires a shapefile(polygons with values) and .xml file which includes every values for polygons in each bands.

My workflow was;

Creating Polygons with mean & sd values for each band I used in. which was numbered over 30k.
Then I've exported that shapefile into .xml like OGR plugin required from me.

Problem started in this point;
Whenever I try to compute and right after train a model for final prediction, plugin doesn't produce any output even if I've executed it smoothly.

The modules I've used in a row are;

ComputeOGRLayersFeaturesStatistics
TrainOGRLayersClassifier
OGRLayerClassifier

Step 3 doesn't work just because I've got no model from step 2, even if I define the output into a file.

What would be the problem?
I hope I've set the problem clearly.

P.S. there's no warning nor error in the process when I execute it. There's nothing wrong with module Log.
P.SS I'm using QGIS 2.14.7 with all the modules in it.


texts:
I have worked a piece of academic work using the PostGIS extension on PostgreSQL (pgadmin). I have done a lot of work include buffer generation and data processing.

How do I cite this software for the publication purpose?

For example the r software R (R Core Team, 2015)

R Core Team, 2015. R: a language and environment for statistical computing. In: Team
RDC (Ed.), R Foundation for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.


texts: Why aren't my basemap tiles showing up via ArcGIS Server

texts:
I have a leaflet js script referencing the ArcGIS server imagery api.  The code worked fine up until about two days ago.  Now when I click on the layer there's a gray background.  I am using my local machine to deploy the map.  My cloudmade maptiles are appearing without error, and there in no error in the console regarding truncated images or anything of the sort for the ArcGIS imagery api.  I am using this api.
http://www.arcgis.com/home/item.html?id=10df2279f9684e4a9f6a7f08febac2a9

var aLayerTwo = L.tileLayer('http://services.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer',{

                attribution: 'Map data &copy; ESRI, Created @ CSUN Geography',
                });


texts: Using iterator for lasd files in ArcGIS ModelBuilder?

texts:
I am trying to use model builder to create a number of DEMs from LiDAR point clouds. To do so, I need to filter the point clouds for ground points only and then convert the points to raster. In the past, I have used the tool "LAS Dataset to Raster" but I am unable to set the filters through this tool. (I would have to go into the properties of each layer individually and set the filter which is not an option for the number of files I have). In searching other questions, I saw that the tool "Make LAS Dataset Layer" allows me to set the filters and combined with "LAS Dataset to Raster" I can set the filters and create the DEM through model builder (see this question and answer: How to specify LAS dataset filters in ArcMap ModelBuilder?).

My lingering issue is that I cannot get any of the iterators to work with lasd files. And when I try to connect the iterator to "Make LAS Dataset Layer", my only option is as a precondition. How do I connect this process to an iterator with lasd files?


texts:
I am trying to make some calculations using PostGIS and PostgreSQL. There are 3 important tables : wbn, wgo and wvb. These tables have about 500.000 objects in them (geometries). The problem is that my calculations are taking way to much time. I think I found the query's responsible for this;

INSERT INTO bio_schema.wgotemp1(wbn_uidn,wgo_uidn)
            (
                SELECT gdi_agiv_schema.wbn.uidn, gdi_agiv_schema.wgo.uidn
                FROM gdi_agiv_schema.wbn, gdi_agiv_schema.wgo
                WHERE ST_Intersects(gdi_agiv_schema.wgo.the_geom, gdi_agiv_schema.wbn.the_geom)='t'
                    AND ST_Length(ST_Intersection(gdi_agiv_schema.wgo.the_geom, gdi_agiv_schema.wbn.the_geom)) > 2
                    AND gdi_agiv_schema.wbn.uidn = f_uidn
                    AND gdi_agiv_schema.wgo.type = f_type

            );

AND

INSERT INTO bio_schema.tempwvb(wbn_uidn,wvb_uidn,sum)
            (
                SELECT gdi_agiv_schema.wbn.uidn, gdi_agiv_schema.wvb.uidn, ST_Length(ST_Intersection(gdi_agiv_schema.wbn.the_geom, gdi_agiv_schema.wvb.the_geom))
                FROM gdi_agiv_schema.wbn, gdi_agiv_schema.wvb
                WHERE ST_Intersects(gdi_agiv_schema.wbn.the_geom, gdi_agiv_schema.wvb.the_geom)='t'
                    AND gdi_agiv_schema.wbn.uidn = f_uidn
            );

What they are doing is looking for geometries (from the tables wgo and wvb) located in a certain polygon (from table wgo). This takes a lot of time, but for my calculations it is absolutely necessary to know which objects from the tables wgo and wvb are located inside this polygon and what their length is.

What can I do to speed up my queries?

Thanks


texts: How to update a base map to get a better image

texts:
I want to use a base map under some layers I'm working on. It's in an island in Azores, Portugal. I need one World Imagery so It can be seen all the vegetation, buildins and streets. My problem is that those maps have a big flaw right in the place I need. You can see it in the picture. How can I fix it?

I'm using ArcGIS 10.1.


texts: Creating multiple Feature Datasets from table/list using ArcPy?

texts:
How to Create multiple Feature Datasets from a table/list. I can create a single feature dataset in a geodatabase with this Python code.

import arcpy
from arcpy import env

env.workspace = "I:/python/MultipleFD2GDB"
arcpy.CreateFileGDB_management("I:/python/MultipleFD2GDB", "HabitatAnalysis.gdb")
arcpy.CreateFeatureDataset_management("I:/python/MultipleFD2GDB/HabitatAnalysis.gdb", "Cultural_Resources", "I:/python/MultipleFD2GDB/2229.prj")

These are the list of Feature dataset I want in the geodatabase.
"EMI_EMF", "Cultural_Resources", "Parcels", "Hazardous_Materials", "Footprint", "Checkpoint_B", "Wetlands", "Botany", "Land_Use", "Buffers", "Air_Quality", "Transportation", "Topography", "Aesthetics_VisualQuality", "Socioeconomics", "Noise_and_Vibe", "Safety_and_Security", "Alignments", "Annotation", "Mapbook_Grids", "Biology", "Geology", "Public_Utilities", "Agriculture", "Public_Lands", "Hydrology", "BSA", "Overview", "Cumulative", "Alternatives_Analysis", "Permission_to_Enter", "Public_Affairs", "Right_of_Way", "Wells", "Wind_Energy"


texts: Extracting pixel values and their n neighbours from a raster given a set of coordinates

texts:
I recently saw Get surrounding pixel value from point shapefile in R or Python on GIS SE. Although it was a rather broad question and got closed, I found it quite useful and interesting. Thus, I decided to reformulate it and give two possible solutions so future readers can benefit from them.

Given a set of coordinates (e.g. a CSV file, a GeoJSON file or a Shapefile) and a raster (both sharing the same spatial reference, of course), how can one extract the corresponding pixel values and their n neighbours (e.g. 8, 15, 24 or 35) using Python?


texts: Invalid syntax when trying to define a dictionary in pyQGIS code

texts: I'm trying to run a code in which I need to use a dictionary, but the terminal acuses 'Invalid Syntax' when I try to define that dictionary. Outside of QGIS the code runs without any problems.


texts: Getting ST_Difference between two different tables in postgis

texts:
One table has a series of grid squares in it, and the other has millions of voronoi polygons in it. I need to create a table that has the difference between the two of the them so that I effectively get the holes in the voronoi layer.

I've tried several things like this:

SELECT st_astext(ST_Difference(density_squares.geom,
       st_collect(density_squares_voronoi.geom))) As new_geom
FROM density_squares_voronoi left JOIN density_squares ON
       ST_Intersects(density_squares.geom,density_squares_voronoi.geom)
group by
       density_squares.geom;

But I get one output gridsquare for every voronoi polygon instead of the one cropped polygon like I need. Will I need to union the voronoi layer first or something?

I had to crack on and get it finished, so I just wrote some .net code to run it as a batch process. However I'm now trying to union the data and that is also taking a long time. A screen shot is here of what it looks like

and the sql i'm trying to use is this:

Select st_union(density_square_islands.islands)
into density_square_islands_union
from density_square_islands
where st_touches(density_square_islands.islands,density_square_islands.islands)

You can see the individual squares within the cropped islands. I just need to union the squares within each island into a single polygons


texts: Creating series of dialog boxes for parameter variables in sub-models of ModelBuilder?

texts:
I'm working on a model to automatically generate ideal routes for mail delivery. The first thing it needs to do is ask the user whether he/she is doing a full route or a partial route. Then, it should go to one of two separate processes to either generate directions for a full route (based on the user putting in a single route number) or to generate directions for specific coded stops (entered individually by the user).

I have an if/then question at the beginning (Full Run or No) which asks the user whether they are doing a full route (Y or N). However, when I run this as one big model, the initial dialogue box includes the Full Run question (Y or N), a space for the user to enter the full route query, AND the a space for them to enter specific stops.

When I make full run and partial run sub-models, neither shows up in the dialogue box at all and the model runs but doesn't really do anything.

What I would like is a series of sequential dialogue boxes. The first one should ask "Full Run or Partial Run?" and the next one should EITHER be the run number request OR the multivalue stops request (based on whether the user picked Y or N). Anyone know of a way to do this? I know very little Python, so a pure ModelBuilder solution or just a conceptual explanation would be ideal

This is the main model:

And this is what each of the two sub-models look like:


texts: Using an expression on the Label Manager tools

texts:
I am trying to display two labels in each polygon and I want to show One label showing  Year 1 and show another line showing in 140,000 acres. I will use that in the Label Manager and use it in an expression.

I was able to write like this Math.round([Acres]) and it only shows 14000 but does not show the comma like this 14,000.

How do I insert a comma to write in an expression ?


texts: Getting PNG image from geoserver with centered feature

texts:
I want get PNG image from geoserver, there will be centered feature.
Problem is bounding box.

How get i extent one feature from wms layer? I use cql_filter for select one feature.

One solution is calculate bbox from geometry, scale and sizes image (point layer).

There is a better result, some like "getExtent" parameter?


texts: What's the best way to run a visibility analysis in PostGIS 2?

texts:
I am looking for the best way to run a visibility analysis in PostGis 2, i.e. given 2 points and obstacle polygons, are these 2 points visible from each other ?

I did not find anything relevant for now. Are you aware of :

such a function
a library I could use to perform this
or do I have to implement it by myself ?

Thanks in advance


texts: Rounding numbers in Raster Legend in QGIS?

texts:
Is it possible to have the numbers in the raster legend automatically rounded off?

In the example below I'd like the labels to be 8, 29, 51 and so on.

I know I can manually edit the labels, but is it really necessary to do that every time I make a map?


texts:
I have a group of polygons. My objective is to find the duplicate polygons from these group of polygons.

Using NetTopologySuite(A c# port of JTS), it is possible to compare two geometires and check if they are equal. The brute-force method(checking each polygon against every other polygon) is the only idea that comes to my mind but it is not usable if there are large number of polygons. Are there any algorithm's which would better the brute-force in this case?

I found this script  which seems to contain ideas to what i am looking for(the comments indicate a divide and conquer approach..but not much details in the comments). But I have to admit i cannot get anything out of it :) .. It was done for ArcView, which I am not familiar with.

Note: I am not looking for a PostGIS/database solution as in this question. I am looking for something that can be integrated to AutoCAD, Quantum GIS or such desktop GIS products through customization(c#,c++, python etc..)


texts: After making transparency of an orthophoto, I can see thin black square lines in the pdf of a layout? Anyone know if you can fix it?


texts: In QGIS are there similare construction tools to ArcGIS?

texts: For a project we are using ArcGIS to reshape autogenerated features. For that we often use the construction tools polygon, rectangel, circle, etc in ArcGIS. The right angle is also an often used tool. Are there some similar construction tools in QGIS? I found the CADTool-plugin and the Freehand-plugin, but that seems not to be the best solution.


texts: Looping FeatureCollection and ImageCollection in Google Earth Engine editor

texts:
I have an image collection and feature collection. In my image collection, I have filtered the dates and a specific band by this:
'''
var imagecollectionL8 = ee.ImageCollection("MODIS/006/MOD11A1").
                        filterDate('2020-02-01', '2020-06-30').
                        select(['LST_Day_1km']);    //About 150 Images
'''

My featureCollection contains list of counties in United States:
var counitesList= ee.FeatureCollection("TIGER/2018/Counties");

I wanted to loop my imageCollection through this featureCollection and extract "LST_Day_1km" for each county on all my filtered dates and export it to a drive as CSV. I followed scripts in Google editor API but unable to find the correct method(tried clipping, reduceRegions).
I added this function at the end of the code
Export.table.toDrive({
  collection: finalCollection,
  description: 'l8_finalCollection',
  folder: 'finalCollection_test',
  fileFormat: 'CSV'});

Still dint work.
Another approach which I tried is, separating the firs Image in my image collection and passing it through the below code;
var Features = firstImage.reduceRegions   //Just FirstImage
  ({
  collection: counties,
  reducer: ee.Reducer.mean(),
  scale: 30,
  });

  Export.table.toDrive
  ({
   collection: Features,
   description : "vectorsToDriveExample",
  fileFormat: "CSV"
  })

I got the LST_Day_1km and is it the right approach to do it? Or any other better ways?
I want to iterate through my image collection rather than reducing it to a single image. How can I do this with Python API or JavaScript inputs?


texts: Why some of my vector layers do not display and some of them do?

texts: I was trying to re-add the water, building, and protected areas layer again but it does not work, It just displayed landuse area and roads, do you know why this could happen?


texts: How to unit test ArcObjects with Mocking?

texts:
I'm a huge fan of unit testing, but still make use of a FGDB to grab features for running unit tests against when using the ArcObjects framework.

Is anyone successfully using mocking against thinks like IFeature, IGeometry, IWorkspace etc.  If so, I'd love to see some examples of how you're doing it.  I don't really care about what mocking framework you use, just seeing how you're doing it would be greatly appreciated.

The problem I see is that you're having to slice and dice between so many interfaces on the same object, that the overhead of creating a representative mock object would be huge.


texts: Exporting Time Series Data using Google Earth Engine and JavaScript?

texts:
I'm new to the GEE platform.

How can I create a time series analysis (vegetation indices) of a particular shape?

would like to export in table form

 var OLI = ('LANDSAT/LC08/C01/T1_RT_TOA');

//lOCALIZAÇÃO

var shp = ee.FeatureCollection(table);

// DEFINIÇÃO DO CENTRO DA ÁREA DE ESTUDO

Map.centerObject(shp);

// SELEÇÃO IMAGEM LANDSAT EM FUNÇÃO DO TEMPO E LOCALIZAÇÃO

var imagery = ee.ImageCollection(OLI)

var collection = ee.ImageCollection(imagery
    .filterMetadata('CLOUD_COVER','less_than',3)
    .filterDate('2015-01-01','2019-09-17')
    .filterBounds(shp));

// ADICIONANDO AS IMAGENS CLIPADAS NO CONSOLE DIREITO

print ("Coleção de Imagens clipadas: ", collection);

// MOSAICO E CLIP DA IMAGEM

var clipimagem = collection.mosaic().clip(shp);

// ADICIONANDO AS IMAGENS CLIPADAS NO CONSOLE DIREITO

print("Imagens Cor verdadeira: ", clipimagem);

//------------------ CALCULOS ÍNDICES DE VEGETAÇÃO ---------------------------

  // Use the normalizedDifference(A, B) to compute (A - B) / (A + B)
var ndvi = clipimagem.normalizedDifference(['B5', 'B4']);

// Use the NBR (A, B) to compute (A - B) / (A + B)
var nbr = clipimagem.normalizedDifference(['B5', 'B7']);

var evi = clipimagem.expression(
    '2.5 * (nir - red) / (nir + 6 * red - 7.5 * blue + 1)',
    {
        red: clipimagem.select('B4'),    // 620-670nm, RED
        nir: clipimagem.select('B5'),    // 841-876nm, NIR
        blue: clipimagem.select('B2')    // 459-479nm, BLUE
    });

var savi = clipimagem.expression(
    '2.0 * (nir - red) / (nir + red + 1.0)',
    {
        red: clipimagem.select('B4'),    // 620-670nm, RED
        nir: clipimagem.select('B5'),    // 841-876nm, NIR
    });

// ADICIONANDO LAYER AO MAPA

// Make a palette: a list of hex strings.
var palette = ['FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',
               '74A901', '66A000', '529400', '3E8601', '207401', '056201',
               '004C00', '023B01', '012E01', '011D01', '011301'];

Map.addLayer(shp.draw({color: 'ffffff', strokeWidth: 5}), {},'Parque Estadual do Cantão');
Map.addLayer(clipimagem, {bands:'B5,B4,B3'}, 'FALSA COR');
Map.addLayer(ndvi, {min:0, max:1, palette: palette}, 'NDVI');
Map.addLayer(nbr, {min:0, max:1, palette: palette}, 'NBR');
Map.addLayer(evi, {min:0, max:1, palette: palette}, 'EVI');
Map.addLayer(savi, {min:0, max:1, palette: palette}, 'SAVI');


texts:
I calculated the azimuth angle between two points. I'm trying to create a new field where I will store my angle values, but although the field is created it stays empty, no values within!

I'm using ArcGIS 10.1, the Python window. This is the code:

import arcpy
import math
arcpy.AddField_management("Parcelles_class_FeatureToPoi3","new_field","DOUBLE")
cur = arcpy.UpdateCursor("Parcelles_class_FeatureToPoi3")
for row in cur :
row.getValue("POINT_X")
row.getValue("POINT_Y")
row.getValue("end_x")
row.getValue("end_y")
row.setValue("new_field", math.atan2((end_y-POINT_Y), (end_x-POINT_X)))
cur.updateRow(row)


texts: Custom basemap not showing surroundings in CartoDB

texts:
I have entered the following png into the XYZ URL for the Carto custom map but it is not working (for either my link below or their example link). I keep getting my custom map, which I rectified on Map Warper, just floating in blank space:

My rectified map:

The image I see on Carto:

The Carto blog on the topic.

Any ideas why this is happening?


texts: Use arcpy.mapping to list broken data layers?

texts:
I want to list all the broken links in all mxds in a server drive. But whenever I try I get an error come up. Can someone use a specific example, say the drive is D drive, how would I put it into ESRIs script here:

import arcpy, os
path = r"C:\Project"
for fileName in os.listdir(path):
fullPath = os.path.join(path, fileName)
if os.path.isfile(fullPath):
    basename, extension = os.path.splitext(fullPath)
    if extension == ".mxd":
        mxd = arcpy.mapping.MapDocument(fullPath)
        print "MXD: " + fileName
        brknList = arcpy.mapping.ListBrokenDataSources(mxd)
        for brknItem in brknList:
            print "\t" + brknItem.name
del mxd

Also I'm not used to not having pythonwin, does anyone know where to d/l it? I googled it but it kept telling me to install it from the ArcMap cd (I'm just using a work computer).


texts: MapInfo Pro 12.5 halo font in layout window

texts:
I wonder if it is possible to add halo font in Layout window in MapInfo Pro 12.5.

I was trying to do this but this option is disable.
Is any way to do this in MapInfo Pro or how to do this?


texts:
I'd like to reuse a page layout across multiple pages - specifically the page border, the title, etc. Is there an easy way to do this?

Right now I'm making new pages and copy/pasting everything.

Using QGIS 3.4


texts: How to detect/select/identify island polygons

texts:
I am generalizing geological maps and as a part of the process I am dealing with "Island Polygons". For that, I first need to identify them or select them.

I already tried Select by Location, by creating another layer from the same layer and used intersect the source layer, which select not only island polygons but polygons which touch the polygon also.

How do I select only Island polygons which are surrounded only by one polygon?


texts: Getting value coded of domain from map layer with arcgis javascript api

texts: I'm trying to read domain values on a map, how to achieve this objective in arcgis javascript api? any help?


texts: Planet API KeyError:location problem with get_assets_by_id()

texts:
I'm trying to download a list of planet scenes with get_assets_by_id. Scenes were selected visially through Planet Explorer and theirs IDs  were copied.  But I keep recieving messege KeyError: 'location' after the first image in the list have benn dowloaded.

Here is my code:

for ids in id_list:
    my_asset = client.get_assets_by_id('PSScene4Band', ids).get()
    client.activate(my_asset['analytic'])
    my_asset = client.get_assets_by_id('PSScene4Band', ids).get()
    callback = api.write_to_file('**dir**')
    body = client.download(my_asset['analytic'], callback=callback)
    body.await()
    time.sleep(30)

So, what is the reason for this error?


texts: PostGIS: Storing large gridded time series

texts:
I'm trying to store a large raster timeseries on PostGIS-2.2 / PostgreSQL 9.5.  The data comes as a collection of HDF-EOS (basically HDF4) files.  Each file contains 24 bands, representing hourly values of the whole grid for one day.  There is 16 years of data (5844 files).

At present, I'm storing each day as a raster with 24 bands, basically because it's the easiest to implement the import for it.  The raster is split into 12 tiles, 42x45 (this was decided by raster2pgsql's -t auto option).  I then add an extra field to the table for the date of the raster and after importing each file, UPDATE swgdn SET raster_date='XXXXXX' WHERE raster_date IS NULL.

It is, in a word, slow.  So far the import's been running for around 2.5 hours and it's processed nearly 3% of the import data.  I'm not measuring it, but it feels like it's getting slower as more data is inserted.

Is this a reasonable strategy for storing this type of data?  Is there one that will perform better?

Would storing the whole timeseries as bands of a single raster perform better?  That would be 140256 bands (and likely to grow).

For reference, here's a sample of the script I'm using to import the data:

raster2pgsql -I -t auto -c 'HDF4_EOS:EOS_GRID:"file1.hdf":EOSGRID:SWGDN' swgdn | psql
psql -c "alter table swgdn add raster_date text;"
psql -c "update swgdn set raster_date='20000101';"
raster2pgsql -I -t auto -a 'HDF4_EOS:EOS_GRID:"file2.hdf":EOSGRID:SWGDN' swgdn | psql
psql -c "update swgdn set raster_date='20000102';"

and so on.


texts: Python code making thousands of empty features in QGIS?

texts:
Am a used to ArcPy and using that for my coding, however, I now have to use QGIS. So I want to be able to write a script that will take a polygon that I made and then duplicate said polygon two meters over and then do that however many times I want.

Here is an example of an end product would be

So the problem is that this code creates thousand of empty feature but also makes the polygons that I want it to make. How would you solve this?

from qgis.core import *

layer = QgsVectorLayer(r"C:\Users\JBlankenship\Documents\Wilbur-Ellis\Drone Pics 2018\Iowa\7004-1824\June_5\Shapefile\7004_1824.shp", "7004_1824.shp", "ogr")

The_Features = layer.getFeatures()
times_x=0
times_y=0
num1=0
num2=0
num3=0

while num1 <= 7:
    print("Big loop",num1)
    num2=0
    num3=0
    x1=0
    x2=15.4+times_x

    y1=4.6
    y2=0

    for i in layer.getFeatures():
        geom = i.geometry()
        Polygon_cor = geom.asMultiPolygon()
        Point_1_x = Polygon_cor[0][0][0][0]
        Point_1_y = Polygon_cor[0][0][0][1]

        Point_2_x = Polygon_cor[0][0][1][0]
        Point_2_y = Polygon_cor[0][0][1][1]

        Point_3_x = Polygon_cor[0][0][2][0]
        Point_3_y = Polygon_cor[0][0][2][1]

        Point_4_x = Polygon_cor[0][0][3][0]
        Point_4_y = Polygon_cor[0][0][3][1]
        print("Reset")

        break
    print(num1)
    if num1 > 0:
        while num3 < 1:
            num3+=1
            Point_1_x = Point_1_x + x2
            Point_1_y = Point_1_y + y2
            Point_2_x = Point_2_x + x2
            Point_2_y = Point_2_y + y2
            Point_3_x = Point_3_x + x2
            Point_3_y = Point_3_y + y2
            Point_4_x = Point_4_x + x2
            Point_4_y = Point_4_y + y2
            print("Done R",num3)
            times_x += 15.4

            # Set the provider to accept the data source
            prov = layer.dataProvider()
            points = [QgsPointXY(Point_1_x,Point_1_y), QgsPointXY(Point_2_x,Point_2_y), QgsPointXY(Point_3_x,Point_3_y), QgsPointXY(Point_4_x,Point_4_y)]

            # Add a new feature and assign the geometry
            feat = QgsFeature()
            feat.setGeometry(QgsGeometry.fromPolygonXY( [points]))
            prov.addFeatures([feat])

            # Update extent of the layer
            layer.updateExtents()

    while num2 < 7:
        num2+=1
        Point_1_x = Point_1_x + x1
        Point_1_y = Point_1_y + y1
        Point_2_x = Point_2_x + x1
        Point_2_y = Point_2_y + y1
        Point_3_x = Point_3_x + x1
        Point_3_y = Point_3_y + y1
        Point_4_x = Point_4_x + x1
        Point_4_y = Point_4_y + y1
        print("Done C",num2)

        times_y += 1

        # Set the provider to accept the data source
        prov = layer.dataProvider()
        points = [QgsPointXY(Point_1_x,Point_1_y), QgsPointXY(Point_2_x,Point_2_y), QgsPointXY(Point_3_x,Point_3_y), QgsPointXY(Point_4_x,Point_4_y)]

        # Add a new feature and assign the geometry
        feat = QgsFeature()
        feat.setGeometry(QgsGeometry.fromPolygonXY( [points]))
        prov.addFeatures([feat])

        # Update extent of the layer
        layer.updateExtents()

    num1+=1


texts: Converting stream definition threshold to extract drainage network to metric scale?

texts:
I extracted two drainage networks from a DEM of 30 meters resolution with two stream definition thresholds: 1000 pixels and 500 pixels.

Can you tell me what this thresholds represent in terms of scale (1/25000, 1/50000...etc)?


texts: Should 2D Positions and 3D Positions be considered conceptually different entities?

texts:
Let's suppose we want to model GeographicPositions, described by coordinates or coordinate tuples inside some Coordinate Reference System. Also, these would be put inside some sequence to represent Lines (Polylines, Tracklogs, etc.)

If we consider a Mapping Software Application, where we can draw over a map image with the mouse (such as Google Maps, OziExplorer, etc.), we are defining horizontal positions based on Latitude and Longitude, so these positions are inherently 2D.

Now if we consider a GPS device, it tracks spatial positions based on Latitude, Longitude and Elevation, so these positions are inherently 3D.

We could, of course, use DEM to get terrain elevation based on horizontal position, but that should not (as I understand) be considered an INTRINSIC PROPERTY of that position, since it is defined by an external function (an operation over the coordinates AND the DEM dataset).

In the other hand, a GPS device can be moved not only along the terrain surface, but it can be inside an airplane, and thus its elevation is intrinsically defined and independent of horizontal position, not being "extractable" by any operation over DEM dataset.

So, my questions are:

Is it common practice in GIS Systems to treat 2D and 3D positions (and collections of positions) as conceptually different entities?
Should a 2D position be considered a limited type of a more general 3D position?

With Elevation set to None/Null?
With Elevation set to 'NaN'?
Without Elevation property at all?

Should a 3D position be considered an extended type, derived from a more basic 2D position by adding an Elevation property?
How to represent a 3D tracklog where some of its positions' elevations are missing? How should it be mathematically processed regarding its missing values? Should a tracklog instance like this exist in the first place?

I would be glad to update the question if anyone thinks it must be clarified!


texts:
I currently am developing a map that uses the following code to set the map background:

//set the map background
L.tileLayer('http://{s}.tile.thunderforest.com/transport/{z}/{x}/{y}.png', {
    attribution: 'http://www.thunderforest.com/'
}).addTo(map);

I have signed up for a Nokia HERE account so that I can use their transit map. I have an APP ID and APP CODE.

According to their instructions, (which I have also included below), I need to do a few steps before I can use their maps.

Can someone help me translate their steps into the what I need to do?

// Create a Platform object (one per application):
var platform = new H.service.Platform({ 'app_id': '{YOUR_APP_ID}',
    'app_code': '{YOUR_APP_CODE}' });

// Get an object containing the default map layers:
var defaultLayers = platform.createDefaultLayers();

// Instantiate the map using the normal map as the base layer:
var map = new H.Map(document.getElementById('mapContainer',
  defaultLayers.normal.map);

// Change the map base layer to the satellite map with traffic information:
map.setBaseLayer(defaultLayers.satellite.traffic);


texts:
I am new to QGIS and GIS in general. I have a CSV file that looks like this:

FIPS Code   Area Name   Area Description
10001       A1          D1
10002       A1          D1
10003       A2          D2
10004       A2          D2

I want to create a vector layer or shapefile that will show boundaries for each Area based of FIPS codes that belong to that Area. I do not want to see county FIPS borders, just Area borders.

I have "wkt_geom" data for each FIPS. I extracted it from Counties Vector Layer from ESRI's website. It looks like this:

MultiPolygon (((-99.4097266969388329 28.63827052216430502, -100.11285992001091927 28.63845898525977418, -100.11502942137985883 29.08920919633010271, -99.41330138287827367 29.09272057773335973, -99.4097266969388329 28.63827052216430502)))

Not sure if this is helpful.


texts: Get coordinates from drawn items with Mapbox-gl-draw

texts:
I'm using mapbox-gl-draw in my application and I need to get and save the coordinates of the drawn items on the map. I've checked and I found turf getCoor method, but I don't know how to access the info because I'm using the default set of tools, which you obtain thru the lines below

var Draw = mapboxgl.Draw();

map.addControl(Draw);

Can someone help me figure this out? I'm sure it's pretty easy, I just can't find as many examples online as I wanted it.

EDIT

Also, while trying to figure this out I got type error from turf, which it means just calling <script src='//api.tiles.mapbox.com/mapbox.js/plugins/turf/v1.3.0/turf.min.js'></script> in my html file didn't work. Am I including it wrong? with https: before the slashes didn't work either.


texts: What does it mean when ST_Equals returns false and `=` returns true?

texts:
There are two methods of determining equality in PostGIS.

The equals operator x.geom = y.geom
ST_Equals(x.geom,y.geom)

Obviously these two methods work different, but the docs on ST_Equals are not clear why,

  Returns TRUE if the given Geometries are "spatially equal". Use this for a 'better' answer than '='. Note by spatially equal we mean ST_Within(A,B) = true and ST_Within(B,A) = true and also mean ordering of points can be different but represent the same geometry structure. To verify the order of points is consistent, use ST_OrderingEquals (it must be noted ST_OrderingEquals is a little more stringent than simply verifying order of points are the same).

This seems to mean there is two criterion for ST_Equals, and that is that

ST_Within(A,B) = true
ST_Within(B,A) = true

However, that ordering of points does not matter with ST_Equals. Presumably, that would matter with the equal operator then =? But, that doesn't match up with my observation.

So as three in the same questions..

What does the equal operator (=) do?
What does the ST_Equals() do?
What criteria does spatially equal include that = does not include?


texts: openlayers python proxy.cgi script is not taking parameters

texts:
I have setup the python proxy.cgi in cgi-bin folder in my server.

[myserver]/cgi-bin/proxy.cgi returns the openlayers website

[myserver]/cgi-bin/proxy.cgi?url=www.openlayers.org returns list index out of range (www.openlayer.org is defeined in the allowed host)

[myserver]/cgi-bin/proxy.cgi?url="www.openlayers.org" returns list index out of range (www.openlayer.org is defeined in the allowed host)

[myserver]/cgi-bin/proxy.cgi?url=http://www.openlayers.org returns page not found

[myserver]/cgi-bin/proxy.cgi?url="http://www.openlayers.org" returns page not found

The proxy is working, but if i use url without http it gives list index out of range, if i use http in url the script is failing.

The url i wasnt to access is http://www.openstreetmap.org/api/0.6/way/190554764 and the cross domain issue is not allowing.


texts: Data providers equivalent to USDA/NRCS/FSA Common Land Unit (CLU)?

texts:
For one of my agricultural projects, I really need reasonably accurate field boundaries for fields in the US, specifically Iowa right now but later, other similar states.

The USDA/NRSC/FSA Common Land Unit (CLU) data is exactly what I need, but apparently since 2008 it cannot be accessed by the public.

Does anyone know of a similar source, commercial or otherwise, of this type of data? I don't need much for attributes, just field boundaries that will let me aggregate analytic data within a field, accurately.


texts:
There is a method that accepts JSON data. Based on this data, points are created and added to the graphics layer. How to implement search in a graphics layer? I found an example of how to implement search in FeatureLayer (sample)! How to do the same in GraphicsLayer?

A method that draws the JSON points:

function draw_point_layer(stations, list_name_net) {
                console.log(stations);
                for (i = 0; i < stations.length; i++) {
                    var point = {
                        type: "point", // autocasts as new Point()
                        longitude: stations[i]['longitude'],
                        latitude: stations[i]['latitude']
                    };
                    var pointAtt = {
                        Code: stations[i]['sitecode'],
                        Name: "" + stations[i]['name'] + "",
                        Type: stations[i]['type'],
                        Nets: list_name_net[i]
                    };
                    var pointGraphic = new Graphic({
                        geometry: point,
                        symbol: markerSymbol,
                        attributes: pointAtt,
                        popupTemplate: {
                            content: [{
                                type: "fields",
                                fieldInfos: [
                                    {
                                        fieldName: "Code",
                                        label: 'Код',
                                    },
                                    {
                                        fieldName: "Name",
                                        label: 'Имя',
                                    },
                                    {
                                        fieldName: "Type",
                                        label: 'Тип',
                                    },
                                    {
                                        fieldName: "Nets",
                                        label: 'Сети',
                                    },
                                ]
                            }]
                        }
                    });
                    pointLayer.add(pointGraphic);
                }
            }


texts: Unable to choose desired layer when georeferencing in ArcMap?

texts:
I am georeferencing some images in ArcMap 10.3.1.

When I want to add the specific layer I'm going to georeference, I cannot choose this from the drop-down menu.

Any ideas on how to solve this problem?


texts: Cached vector tiles with dynamic attributes

texts:
I am looking to implement a web map of property boundaries and owner information. The parcel data (property boundaries) is updated infrequently, but the assessment data is updated weekly, as properties are bought and sold. What are the best practices to cache a somewhat dynamic, large data set?

Here's a list of some of the features I'd like to show on the map:

Shape (infrequently changes)
Pin (infrequent)
Physical address (infrequent)
Current owner (frequent)
Last sale information (frequent)

I am wondering what the best solution would be to caching the vector features:

Cached vector features, with key and physical location information as the only attributes. Upon interaction in the web map, make a second call to load the current owner/sales information. Would need to regenerate the cache when the parcel/boundary data is updated.
Cached vector features, but with all desired information. Would need to regenerate weekly, or create a process to only regenerate tiles where new information was made available.
Another alternative is to not cache, but already seems foolish, rendering the property boundaries from scratch every request seems wasteful, even though some attributes of the features change week by week.

Does anyone have experience with vector tiling to help me decide? Anyone have good experience or examples of using a static tile map and loading attributes using a second call? I'm using a Linux/PostGIS server to host this application.


texts: NoData specification via IEEE NaN in GeoTIFF for ArcGIS

texts:
I am a programmer exporting 32-bit floating point GeoTIFF rasters from my own application that my users need to read into ArcGIS. Unfortunately I am not an Arc user myself and have zero expertise with it, so please target comments accordingly.

Arc is reading and locating my rasters properly, but I am having problems making it understand which grid cells are empty. Since all positive and negative floating point numbers are legal data values in my rasters, I am trying to use IEEE NaNs as my NoData values, which I understood to be an accepted practice. Furthermore, as my rasters are often dimensionally large but sparsely populated, I'm trying to get maximal compression out of the TIFF PackBits algorithm (which really likes to see sequences of identical bytes) by making an optimal choice with regard to my particular NaN value. Since any 32-bit float whose exponent bits are all 1 is considered a NaN, I'm using 0xffffffff (32 1-bits) as my NaN value since it's both a valid IEEE NaN and is also 4 identical bytes.

Arc seems to treat my 0xffffffff NaN values as low but legal values, i.e., it renders them as black similarly to the actual low legal data values in my file. How can I get Arc to treat this particular constant (which is of course not a value that can be expressed as a floating point number since it's an actual IEEE NaN) as my NoData value and not render it at all?

Not being an Arc user it's quite possible that I'm just not setting the right display/import options within Arc to get what I want. If someone could explain how to do this that may well solve my problem.


texts:
I'm trying to perform regression kriging using the code

creating variogram

m <- vgm(var(soil$Ba), "Sph", sqrt(areaSpatialGrid(slope))/4, nugget=0)
v <- fit.variogram(variogram(Ba~RASTERVALU, data= soil),model = m )
plot(variogram(Ba~RASTERVALU, data= soil),v, main="fitted by lt")

running kriging

slopeuk <- krige(soil$Ba~soil$RASTERVALU,locations= soil, newdata= slope, model=m)

However after the final line of code I get the error

Error in gstat.formula.predict(d$formula, newdata, na.action = na.action,  :
  NROW(locs) != NROW(X): this should not occur
In addition: Warning messages:
1: 'newdata' had 2849126 rows but variables found have 195 rows
2: 'newdata' had 2849126 rows but variables found have 195 rows

Here are some images of the data


texts: Why are 12 hour time values printing as date values with ArcPy?

texts:
I have a field named ORIGDTTIME that contains 12 hour date values such as 06:36 AM.  The field is a date  data type in my feature class, however when I print the values from this field I am not receiving any of the actual values that are in the feature class.  The date field is in a feature class in an Oracle SDE connection. The isue is reproduced in a file geodatabase

Here are the values that are being printed

[datetime.datetime(1899, 12, 30, 13, 50), datetime.datetime(1899, 12, 30, 13, 20), datetime.datetime(1899, 12, 30, 13, 9), datetime.datetime(1899, 12, 30, 13, 45), datetime.datetime(1899, 12, 30, 14, 54), datetime.datetime(1899, 12, 30, 10, 28, 19), datetime.datetime(1899, 12, 30, 12, 48), datetime.datetime(1899, 12, 30, 14, 48), None, datetime.datetime(1899, 12, 30, 10, 36), datetime.datetime(1899, 12, 30, 10, 27), datetime.datetime(1899, 12, 30, 13, 46)].

I am using this script to print the values.

import arcpy
fldName = 'ORIGDTTIME'
fc = 'featureclass'
myList = set([row.getValue(fldName) for row in arcpy.SearchCursor(fc)])
print myList

What should be changed to print the correct values in this script?


texts: Using ListTables when accessing File Geodatabase gives ArcPy IOError?

texts:
I am writing an arcpy script and testing against a file Geodatabase. I can open the Geodatabase workspace and iterate through the "ListFeatrueClasses".

However if I try to "ListFields" I get an IO Error.

I have already tried changing the workspace directory in various ways as suggested by other posts with no success. I also made sure the file Geodatabase is fully accessible and currently its on my local drive.

This is the error:

  Traceback (most recent call last):   File
  "C:\Users\username\Source\Repos\Python-proj1\Python-proj1\openprimegeo.
  py", line 29, in
      fields = arcpy.ListFields(fc1)   File "C:\Program Files (x86)  \ArcGIS\Desktop10.5\ArcPy\arcpy__init__.py", line  1138, in
  ListFields
      return gp.listFields(dataset, wild_card, field_type)   File "C:\Program Files
  (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\geoprocessing_bas e.py", line
  346, in listFields
      self._gp.ListFields(*gp_fixargs(args, True))) IOError: "Test_Polygon" does not exist

  C:\Python27\ArcGIS10.5>

Script:

import arcpy
from arcpy import env

arcpy.env.workspace = "C:/Users/username/Downloads/prime/Prime_Template.gdb/Prime_Template.gdb"

print "starting..."

featureclasses = arcpy.ListFeatureClasses()

for fc in featureclasses:
    print fc

fcList = arcpy.ListFeatureClasses()
print fcList
print "\n\n"
print "List Feature Class fields..."
print "\n\n"
for fc1 in fcList:
    fields = arcpy.ListFields(fc1)
    for field in fields:
        print "{0} is a type of {1} with a length of {2}".format(field.name, field.type, field.length)


texts: Getting lat/lng for latest Wkid:3435, Wkid:102671 and rings

texts:
I've seen similar question as Getting lat/lng from wkid, latestwkid and x y coordinates but not sure how to put rings into ArcGIS server REST endpoint.
<img src="https://i.stack.imgur.com/UbXYb.png" alt="enter    "geometry": {
"rings:"
[
[1197289.1030000001, 1759877.251000002]
[1197292.4299999997, 1759549.5487499982]
[1197005.0549999997, 1759549.939000003]
[1196998.6702499986, 1759968.4439999983]
[1197288.1697500013, 1759969.0909999982]
[1197289.1030000001, 1759877.251000002]
],
"geometryType": "esriGeometryPolygon"} image description here" />
Should I take every two coordinates and process separately? I am not really sure what output input/output spatial reference should be and why. Is there any Python program/package/code that can do it?


texts:
I wish to make an update using CASE, on a subset of a table. This subset is based on a spatial query.
It looks like:

update t1
set f1 = CASE WHEN...
         WHEN...
         END
FROM (select * from t2 join t1 on st_intersects(t2.geom, t1.geom)
WHERE *a_field_from_t2* = *a_value*) as foo ;

This update always works with the whole set of t1 i.e. the spatial is not used.

I can't find the reason why. Any help?


texts:
I'm looking to export all the paths that comprise a long hiking trail in Japan. Unfortunately while some of the paths are named, many of them aren't and some are named inconsistently.
Is there a way to manually select OSM features for export, instead of downloading the area within a large bounding box?


texts:
I have always been curious how other folks are storing their spatial data.

Outside of an enterprise environment where much data is stored in enterprise geodatabases, is it better to store data per project in project folders, (The data may be important for other projects too, or maybe for a stand alone project) OR is it better to file everything in specific geodatabases  or main folders - with shapefiles (accessible for many projects)?

How do you manage data accessability?

What has worked well for you?


texts:
I have recently installed QGIS 2.8.2 on my Windows 8.1 laptop with the stand alone .exe installer.  The processing toolbox is missing the GDAL/OGR, GRASS, QGIS algorithms, Saga and the Orfeo Toolbox from the processing toolbox which I have on my 2.8.2 installation on a Windows 7 machine using the same installer file.  I have checked the plugin options and all of them are set to activated.  Also the paths are set too and they were pointed at C:\Porgra~1\Wein\apps\"module folder"

I have tried changing it to what is recommended in Joseph's comment to the .qgis path under the users directory but with no change to what is shown

Has anyone else encountered this problem?


texts:
I have a CSV file with 800 rows and 20 columns that I would like to join with a shapefile (polygon) that has only 10 rows and 1 column (Id).

How should I go about this?


texts:
In Esri's area calculation sample of IGPFunction implementation they use IGPUtilities but they don't release the instance.

On the singletons page in SDK help they list GPUtilities as a singleton object.  So shouldn't our instances of IGPUtilities be released? And also should we not use the Activator class to create the instance?

In a IGPFunction implemented class do people have any preference on creating a class instance of IGPUtilities (like in the Esri Area Calculator example) or create them in the methods when needed?


texts: ArcGIS Online custom pop-up attributes displayed by a query

texts:
While using ESRIs ArcOnline I am building a web application with multiple custom pop-ups. When I click on a building I want to provide the exact zoning code information with a link to the appropriate website full of zoning codes.

All of the zoning information (along with its appropriate links) has been loaded into fields of the building attribute table. I fully understand how the pop-ups work, but I was wondering if there was anyway possible to change what is displayed on the pop-up through a query.

--> CURRENT pop-up display

--> QUERY purpose: the queries purpose would be to change the pop-up display so if "Historic Area?" is "Yes" then the link is displayed. If the answer is "No" then no link is displayed.

In conclusion, I simply want to make the pop-up cleaner, removing "If yes - more info: " and an unnecessary broken link accompanied if the answer is "No".

So is this at all a possible task?


texts: Fixing problem with processing component

texts:
How can I fix the problem with the processing component?
Since it appears to me that it is broken and I need to activate some functions of the Raster option. When I want to activate the processing plugin, a legend appears saying that this plugin is broken, which may be due to the program update and I would like to know how to fix this problem. I am using QGIS version 3. 10. 13


texts: Computing mean of points using ArcGIS Desktop?

texts:
I have a polygon and inside it are points with z-values. I want to get the mean of the z-values within the polygon.

How can I do that in ArcGIS Desktop 10.5?


texts:
I'm using the following to get a list of dates in a field called "EFF_DATE" from an attribute table for some pre-selected features.

selectDates = [row[0] for row in arcpy.da.SearchCursor(selects, "EFF_DATE")]

I then use max() to get the most current date from that list:

selectDates = max(selectDates)

This is my result for In[]:selectDates

  datetime.datetime(2012, 12, 18, 0, 0)

If I run In[]:print(selectDates), I get

  2012-12-18 00:00:00
  That looks better but what I want is just this
  '12/18/2012'
  How do I convert it?


texts:
I start with Create New Shapefile Layer, I give it a name and then select Polygon.
I then create my Polygon shape, its all good to this point. I now show a polygon with an outline and light red fill. From here I get stuck. I cant change the fill or do anything else with it. I think I'm missing a step.b


texts: GraphicsLayer not redrawing visually while a Graphic is being dragged (JS API 4.7)

texts:
I've implemented a 'route between points' site using the 4.7 JavaScript API. The logic works fine, but I'm having a graphical display issue: when a point is picked up and dragged to a new location it doesn't redraw intermittently, the mouse simply drags an invisible point to the new location, where it reappears.

The 'update' event fires every time the mouse moves, but the graphics layer won't redraw on the screen, even when cloning the graphic, deleting it from the graphics layer, and re-adding it - the usual fix for a graphic not re-drawing upon being edited.

Oddly, everything works fine if the user zooms in or out of the map while gripping the point. Things start drawing, and the greyOutAllRouteLines function (which sets the colour of every line on the map to grey, to indicate that with the shifted point they're now out of date) works as advertised. In both cases the remove/clone/re-add method is required, otherwise nothing redraws regardless of map interaction.

It feels like the graphics layers are being updated quite happily while the 'update' events are firing, but they aren't being prompted to re-draw until the map is interacted with, at which point they function as expected. Interacting with the map using zoom functions and similar don't trigger the effect.

Here's the event handler:


texts: Check if layer is a QGSVectorLayer from Python Plugin within QGIS

texts:
I am trying to extend some concepts in the QGIS Workshop tutorial.  Basically, I would like to check if the layer is of type "QGSVectorLayer".  In the QGIS Python console, I am able to type

canvas = qgis.utils.iface.mapCanvas()
cLayer = canvas.currentLayer()
type(cLayer).__name__

with no problem.  However, when attempting the same from the plugin, I get the following error message:

  if type(self.cLayer).name == "QgsVectorLayer": NameError: global
  name 'cLayer' is not defined

The code is as follows:

def handleLayerChange(self, layer):
        self.cLayer = self.canvas.currentLayer()
        if self.cLayer:
            if type(self.cLayer).__name__ == "QgsVectorLayer":
                self.provider = self.cLayer.dataProvider()


texts: Is it possible to define a custom CRS for a shapefile in QGIS that uses hundredths of a US foot?

texts:
I have several polygon shapefiles whose coordinates are defined as 100x their original coordinates projected from the Washington State South HARN State Plane in us-ft (EPSG 2927). I would imagine the purpose of this scheme is to represent more precise coordinates in integer form. Unfortunately, I am having a hard time representing this custom coordinate system as a proj4 specification in QGIS. Each of these shapefiles has a .prj file and no .qpj file.

Things I have tried or thought of:

creating a custom proj file by editing the +units parameter to "us-ft/100" (not defined).
the .prj file has a parameter called "scale-factor" that I tried to set to 0.01, but QGIS does not appear to read this (not sure if the .prj file is being ignored or just this parameter).
using a geoalgorithm to "rescale" each feature in the shapefile by 0.01, so that they actually fit within the standard EPSG 2927 CRS.
Going back to the original application (Geomedia) to reproject the shapefile, which I am finding is difficult with shapefiles. I would also just like to know if there is a more "open source" way of doing this.

What is the best way to accomplish this task in QGIS?


texts: Creating unique IDs for file geodatabase using ArcGIS for Desktop?

texts:
I use ArcEditor licence level, ArcGIS 10.0 version, and we're talking about a file geodatabase.

I have several feature classes (electric energy domain) and I want to connect them using relationship classes, so I want to use their ID's as foreign keys.

I want to automatically create a significant unique ID for each feature. It'a a text type ID, which would be composed of two letters and several numbers (ex. CH00123, CH00124). The ID would give me three info : what type of feature we're talking about (channel- CH), on which street it is (00), and it's unique code (123). I don't want to type this code for every feature I create, I need something to automatically assign an ID to each feature.

I am not sure how to handle this situation. Do I need to use a data integration software? an engine to generate my special ID's?


texts:
I am using osmnx package to visualize a street network. I would like to visualize both (two) lanes in opposite directions of a road. I am using this code https://github.com/gboeing/osmnx/issues/162
How do I add the output point on the map?
import osmnx as ox

place_name = 'Cergy, France'
G = ox.graph_from_place(place_name, network_type = 'drive')
lines = []
for u, v, data in G.edges(keys=False, data=True):
    if 'geometry' in data:
        # if it has a geometry attribute (a list of line segments), add them
        # to the list of lines to plot
        xs, ys = data['geometry'].xy
        points = list(zip(xs, ys))
        #parallel shift distance
        h = 1

        if not data['oneway']:
            # for each point excluding the start point and end point, shift point based on
            # line to next point
            transformed_points = [points[0]]
            # get pairs of points on the line segment
            for p1,p2 in zip(points[1:], points[2:]):
                (x1,y1) = parallel_point_shift(p1,p2,h)[0]
                transformed_points.append((x1,y1))

                transformed_points.append(points[-1])
                points = transformed_points

            lines.append(list(points))
    else:
    # if it doesn't have a geometry attribute, the edge is a straight
    # line from node to node
        x1 = G.nodes[u]['x']
        y1 = G.nodes[u]['y']
        x2 = G.nodes[v]['x']
        y2 = G.nodes[v]['y']

    if not data['oneway']:
        ((x1,y1), (x2,y2)) = parallel_point_shift((x1,y1),(x2,y2), h)

    line = [(x1, y1), (x2, y2)]
    lines.append(line)


texts: How to get rid of layer created by Make Feature Layer (Data Management)

texts:
I like to run some scripts from mxd. Very often I am using Make Feature Layer tool to facilitate queries.

Tool help says:
“Creates a feature layer from an input feature class or layer file. The layer that is created by the tool is temporary and will not persist after the session ends unless the layer is saved to disk or the map document is saved”.

Question:

How to remove newly created layer without closing mxd?

Note arcpy.mapping.RemoveLayer is not working on layer created that way or maybe I am doing something wrong here


texts: How to get rid of isolated pixels in a raster using gdal, numpy functionalities in python ?

texts: I have a raster which I got after passing it through a certain thresholding to map the debris cover pixels. However, there are certain isolated pixels here and there probably due to limitations in the imagery resolution. How do I get rid of these isolated pixels in a python environment using pure gdal and numpy functionalaties?


texts: Clipping raster layer supplied from web (e.g. OSM)

texts:
When I try (in QGIS 3.8 "Zanzibar") to clip a "web layer", i.e. a raster layer provided from the web (e.g. an OSM map loaded via the QuickMapServices plugin) by the usual method
Raster > Extraction > Clip Raster by mask layer or by extent,
the gdalwarp command as generated by QGIS first reports:

Execution completed in x.yz seconds
Results:
{'OUTPUT': 'C:/Users/nn/Desktop/ClippedLayer.png'}

but in fact fails.

The Log output contains:

INFO    GDAL execution console output
FAILURE: No target filename specified. [even though it *has* been]
Usage: gdalwarp [--help-general] [--formats]
[...]
'url' is not recognized as an internal or external command, operable program or batch file.
'zmax' is not recognized as an internal or external command, operable program or batch file.
'zmin' is not recognized as an internal or external command, operable program or batch file.

I suspect that the problem is with the URL-type string coming from the "web layer" URL:
If the layer to clip exists as a real file, the gdalwarp parameters are e.g.

...
-cutline C:/Users/nn/Documents/QGIS/TEST/areaFinal_UTM.shp
-crop_to_cutline C:/OSGeo4W64/nn/TEST/terrain_UTM.tif
...

and all works well.

For a "web layer" it looks like

...
-cutline C:/Users/nn/Documents/QGIS/TEST/areaFinal_UTM.shp
-crop_to_cutline type=xyz&url=hxxp://tile.openstreetmap.org/%7Bz%7D/%7Bx%7D/%7By%7D.png&zmax=19&zmin=0 ["hxxp" substituted for "http"]
...

which (unsurprisingly) fails.

The only way I found to "pseudo-clip" such a "web layer" was through generating an image from it via
Project > Import/Export > Export Map To Image,
where the extent can then come from another layer (e.g. a mask defining the extent).

However, this export dialog requires input of size and/or resolution, which in turn depend on the scale at which the map is currently viewed in QGIS, so it is not really a good general solution.

Hence the questions:

How can a "web-supplied" raster layer be properly clipped to a limited extent in QGIS?
How would the name of such a "web layer" be properly passed to the crop_to_cutline argument?
(This would require running directly from a console, as QGIS 3.8 does not allow editing of its generated GDAL commands any more, alas.)


texts: Using field variable to calculate values in ModelBuilder?

texts:
I am using ModelBuilder and trying to get a user defined field from a user defined polygon.
From this, I want to later populate a different field with calculated values in a model.

Is there a way to get a list of fields from the user defined polygon and use that value to populate the other field?

How would I go about using the user defined field in a VB expression within the calculate value tool?

newfield = [%field%]?


texts: Classification during DEM generation for Bare Earth Surface-ERDAS

texts:
I am using ERDAS photogammetry tools(eATE) to generate a bare earth type DEM from two Pleiades satellite images(pan-sharpened 0.5 cm resolution). I have run the process and the results were just OK. Quite a few of the trees were still in the DEM, showing up as abnormal elevation values in the DEM(say the surface was 200m, one pixel would be 200 m and one would be 205 due to the tree). The stereo Pleiades images were taken in a semi-arid area with brown grass and roughly 6 meter trees that are spread out.  The trees aren't in every pixel but do create alot of noise in the DEM.

Any advice on how to better remove the trees during this Classification step in eATE? I've played with the "greeness" parameter but it didn't help much.


texts: Logical operators and looping back through a process in Modelbuilder

texts:
I am using a series of models in ArcGIS Pro 2.1 to manage gaps and slivers on a very large dataset, and have hit a wall trying to use a logical operator.

I have split the data into multiple smaller datasets within a File Geodatabase and need to iterate through the FGDB to run an Eliminate on any feature classes that have polygons with an area under 0.1 ha.

The model (snapshot attached) uses the If Selection Exists logical operator to continue to the Eliminate process if a selection expression of <0.1 ha exists however the model stops completely if the output is false (i.e. not all feature classes have polygons under 0.1 ha).

How do I get the model to continue to iterate through all the datasets and Eliminate where required without stopping when it hits a dataset that doesn't have polygons under 0.1ha?


texts:
I have been installing GeoServer locally many times for the past year on Windows XP machines. Now I am trying to install it on our new Windows 7 Machines. However the installation (of the stable version 2.5.2) seems to proceed differently, No longer do you get a menu on the start menu to launch it, in fact you have to hunt around to find GeoServer.

When I tried to launch it (from the startup batch file I found in the Bin directory) nothing happened, various re-installations of all the different versions I could download and I found that putting geoserver in the programs directory (as opposed to the default program files (x86) directory) resulted in a installation that worked, once... but now the batch file no longer works. I have varied the port between the default 8080 and various others to no avail, so it does not seem to be a port conflict. I am really confused as to why the menu is not there, the program seems suddenly very different.

Am I doing something wrong?


texts: Fitting the plotted pixels to the country border in QGIS

texts:
I have created a plot in Python 3.7 over Alaska. Firstly, the pixels covered a lot more than shown in the image (below). Now I have set the lat lon extend to fit the plot 'better'.

Is there a way to optimize this, so the plot actually fits within the border of Alaska?
Can I do it prior to applying it in Python, hence in QGIS?

Here is the Python script
import netCDF4 as nc
import xarray as xr
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature

ds = nc.Dataset(path + 'test.nc')

lons = ds.variables["lon"][:]
lats = ds.variables["lat"][:]
prob = ds.variables["test"][:]

lon_0 = -140
lat_0 = 65

m = Basemap(width = 5000000, height = 3500000, resolution = 'l', projection = 'stere', lat_ts = 40, lat_0 = lat_0, lon_0 = lon_0)

lon, lat = np.meshgrid(lons, lats)
xi, yi = m(lon, lat)

# Plot Data
cs = m.pcolor(xi,yi,np.squeeze(prob))

# Add Grid Lines
m.drawparallels(np.arange(-80., 81., 10.), labels=[1,0,0,0], fontsize=10)
m.drawmeridians(np.arange(-180., 181., 10.), labels=[0,0,0,1], fontsize=10)

# Add Coastlines, States, and Country Boundaries
m.drawcoastlines()
m.drawcountries()


texts: QGIS couldn't load processing plugin - Windows 7

texts:
I'm having a problem loading QGIS.  Usually runs fine on my computer.
Error pops up "Cant'load processing plugin...."  On QGIS 2.16 I was able to get the full error message.  I've fresh installed 2.18.3 and now get this...

"Plugin processing: The plugin will be disabled because it crashed QGIS during last startup. Please report an issue and re-enable the plugin when the problem has been solved."

Any ideas?

The full error is shown below...

Couldn't load plugin processing due to an error when calling its classFactory() method

ImportError: cannot import name uic
Traceback (most recent call last):
  File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python\qgis\utils.py", line 333, in startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing__init__.py", line 37, in classFactory
    from processing.ProcessingPlugin import ProcessingPlugin
  File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python\qgis\utils.py", line 607, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing\ProcessingPlugin.py", line 38, in
    from processing.gui.ProcessingToolbox import ProcessingToolbox
  File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python\qgis\utils.py", line 607, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing\gui\ProcessingToolbox.py", line 31, in
    from qgis.PyQt import uic
ImportError: cannot import name uic

Python version: 2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]
QGIS version: 2.18.3 Las Palmas, 77b8c3d

Python Path:
C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python
C:/Users/msm162/.qgis2/python
C:/Users/msm162/.qgis2/python/plugins
C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\matplotlib-1.3.1-py2.7-win-amd64.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\nose-1.3.3-py2.7.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\tornado-4.0.1-py2.7-win-amd64.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\backports.ssl_match_hostname-3.4.0.2-py2.7.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\certifi-14.05.14-py2.7.egg
C:\PROGRA~1\QGIS2~1.18\bin\python27.zip
C:\PROGRA~1\QGIS2~1.18\apps\Python27\DLLs
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\plat-win
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\lib-tk
C:\PROGRA~1\QGIS2~1.18\bin
C:\PROGRA~1\QGIS2~1.18\apps\Python27
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\PIL
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\jinja2-2.7.2-py2.7.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\markupsafe-0.23-py2.7-win-amd64.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\pytz-2012j-py2.7.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\win32
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\win32\lib
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\Pythonwin
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win-amd64.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\wx-2.8-msw-unicode
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\xlrd-0.9.2-py2.7.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\xlwt-0.7.5-py2.7.egg
C:/Users/msm162/.qgis2//python
C:\Users\msm162.qgis2\python\plugins\mmqgis/forms
.


texts: Why does Arcobjects IEditor operation not fire IEditEvents.OnCreateFeature?

texts:
I built an Arcobjects add-in in the past that records the username and datetime every time a feature is created or modified:

Sub OnStartEditing_Handler()
    AddHandler Events.OnCreateFeature, AddressOf OnCreateChangeFeature_Handler
    AddHandler Events.OnChangeFeature, AddressOf OnCreateChangeFeature_Handler
End Sub

Sub OnCreateChangeFeature_Handler(ByVal obj As ESRI.ArcGIS.Geodatabase.IObject)
    Dim feature As IFeature = CType(obj, IFeature)

    'Update Date field
    'Update Username field

End Sub

I have recently been working on another add-in that updates a different field during an edit session:

'Populate Feature attribute
m_editor.StartOperation()
m_feature.Value(m_hgtFieldIdx) = CInt(height)
m_editor.StopOperation("Measure Height")

The second add-in successfully updates the field, however the first add-in is not updating the username and date fields.  Does the StartOperation() or StopOperation() not fire the OnChangeFeature event?  I cannot find any documentation for this situation.  If not, can I fire it programmatically, or is there another way to catch feature(s) changed during the Start/StopOperation()?


texts: Selecting polygons by raster RGB values using ArcGIS Desktop?

texts:
I'm looking for solutions to this in ArcMap 10.2.2 Basic, with no Spatial Analyst extension.

I am trying to find parcels that have been damaged by a large fire. I have a raster image showing the burn path and surrounding areas. I'm overlaying a parcel layer on top of the raster and would like to extract polygons on top of raster cells that have a specific RGB range.

Here is the image:

The dark red and black pixels always have a red value > 180 and < 210 with green and blue values less than 50. The lighter pixels all have RGB values > 200.

The raster has no attribute table as it is a multi band image. I converted it to a single band but the values are no longer discernible in the burn path from other areas of the image.

I then tried converting the raster to points. Still, no distinct values for the burn path.

Is there a way to select polygons using a specified pixel RGB range?


texts: Control feature rendering order in QGIS with plugin

texts:
I am trying to find the python equivalent of enabling the checkbox for controlling the feature rendering order based on an attribute of the features of a layer (see below screenshot).

I want to base the rendering order on 'height' as shown in the screenshot.

My code below reads:

lyr_features = sorted(lyr_features, key = lambda i: i['height'])
colours = self.linear_gradient(len(lyr_features))
for i, feat in enumerate(lyr_features):
    myMin = feat['height']
    myMax = feat['height'] + 0.5
    myLabel = '{}'.format(feat['height'])
    myColour = QtGui.QColor(str(colours[i]))
    mySymbol = QgsSymbol.defaultSymbol(lyr.geometryType())
    mySymbol.setColor(myColour)
    mySymbol.setOpacity(myOpacity)
    myRange = QgsRendererRange(myMin, myMax, mySymbol, myLabel)
    myRangeList.append(myRange)

myRenderer = QgsGraduatedSymbolRenderer('', myRangeList)
myRenderer.setMode(QgsGraduatedSymbolRenderer.Pretty)
myStyle = QgsStyle().defaultStyle()
defaultColorRampNames = myStyle.colorRampNames()
ramp = myStyle.colorRamp(defaultColorRampNames[25])
myRenderer.setSourceColorRamp(ramp)
myRenderer.setClassAttribute(myTargetField)

lyr.setRenderer(myRenderer)
QgsProject.instance().addMapLayer(lyr)
lyr.triggerRepaint()

This is okay, but I then need to manually go in and enable the feature rendering control. Maybe something like this:

myRenderer.orderByEnabled(True)

This doesn't work but maybe along the lines of this?


texts:
I have a layer like,

realTimeLayer = new esri.layers.GraphicsLayer();

I have an infoTemplate like,

infoTemplate=new esri.InfoTemplate("/*Content goes here*/");

I am creating a graphic, adding it to the layer and plotting the point like,

 var graphic= new esri.Graphic(candidatePoint,candidateSymbol,attr,infoTemplate);
 realTimeLayer.add(graphic);
 map.addLayer(realTimeLayer);
 points.addPoint(candidatePoint);

In this case the info template works fine.
But during some other event, I am adding the graphic like

var graphic= new esri.Graphic(candidatePoint,candidateSymbol,attr,infoTemplate);
map.graphics.add(graphic);
points.addPoint(candidatePoint);

But this time, I don't get the info template in IE9 even after clicking on the push pin. What could be the reason? And how should I correct this?


texts: FindLabel not Defined in ArcMap labeling with Python Parser?

texts:
I'm trying to make an advanced label expression using the python parser and I keep getting the following error:

The expression contains an error.
Modify the expression and try again.

Error 0 on line 0.
Error running expression: FindLabel()
Traceback (most recent call last):
 File "<expression>", line 1, in <module>
NameError: name 'FindLabel' is not defined

def FindLabel ( [Direction] ):
  return  "{}° {}' {}\"".format( [Direction].split("-"))

Can anyone tell me why FindLabel is not defined?


texts: How to get cell value from one raster within a raster catalog using ArcPy?

texts:
I need to be able to obtain the cell value of one cell in one specific raster that is stored withing a raster catalog in a file geodatabase.  The Get Cell Value tool doesn't appear to be able to do this on a Raster Catalog.

Using python, how does one go about doing this?


texts: Shapefile to TAB preserve table structure

texts:
I generally work in ArcGIS and only use QGIS when dealing with TAB files.

I have a client who has very particular attribute table structure requirements and works exclusively in MapInfo.  I thought that I would be able to simply create all of the data in ArcGIS, open the GeoDatabase in QGIS and then save as a .TAB.

However, I have since discovered that when I do this the data and metadata are transferred correctly (with the exception of dates) but the table structure changes - the headings remain correct but their characteristics (mainly length) change to what I am guessing is a default 256.  This has caused some issues...

I have tried using ORG2GUI and ORG2OGR to do the conversion but the same thing appears to be happening.

how do you transform a .GDB or .SHP to a .TAB while preserving the table structure?

I have a very poor work-around using MapInfo 8, lots of copying and pasting, and some finger crossing.  The main issue is that I do not have a more recent version of MapInfo to check that the data is correct before sending and I do not want to pay $6600 for an ArcGIS extension which I think should be part of their basic package.


texts: Failed to Copy SDE Dataset Error 000260 Using arcpy.Copy_management and arcpy.env.overwriteOutput = True

texts:
I thought I could use the arcpy.env.overwriteOutput = True to overwrite the existing datasets within my target SDE "JeffCoGIS.sde"; however, I get an error that states:

ExecuteError: ERROR 000260: Failed to copy Database Connections\GISData.sde\GISData.DBO.GIS_Dept into Database Connections\JeffCoGIS.sde\GISData.DBO.GIS_Dept
ERROR 000260: Item with the same path name already exists.
Failed to execute (Copy).

I am new to python and couldn't find any help in ESRI Resources on how to overwrite the datasets.  Some have mentioned to use replication; however, for our particular workflow, we need to pursue overwrite and copy, if I can get it to work.

Is there a different tool I should use?

If I run this script on the empty SDE, it works fine - just when trying to overwrite the existing datasets.

#Name:  Copy_JeffCoGIS_Rev2.py
#Copy all datasets including the feature classes within each dataset from
#    SDE database "GISData.sde" to SDE database "JeffCoGIS.sde"
import arcpy
arcpy.env.overwriteOutput = True
arcpy.env.workspace = "Database Connections\\GISData.sde"

fcList = arcpy.ListDatasets()
#print fcList

for dataset in fcList:
    arcpy.Copy_management(dataset, "Database Connections\\JeffCoGIS.sde\\"+
    dataset)


texts: Sum up area left "underneath" polygons in QGIS?

texts:
I have an area shape layer and multiple other area shape layers. Now I would like to calculate how much space (square kilometers or percentage of area) of the polygons of this area layer is NOT covered by the other layers. Some of the other layers are overlapping so it isn't as easy as to substract these areas.

Is this possible and if yes - does anyone here know how?

best regards
desputin


texts: How can I get all graphics included into my Feature Layer using ArcGIS JS API 4?

texts:
I'm load my Feature Layer from url and add it to map through property 'layers'. Then I need to highlight polygon and show info popup when click on it. I tried to use hitTest() method in click event from my Scene View, but I get 0 graphics in response. Maybe there is exist another decision to find graphic and highlight it?

Does anybody have any ideas?


texts: How to check hydrological connectivity in flood fill techniques in QGIS?

texts:
I have a set of polygons (include their elevation levels), representing geographical areas. I want to check which polygons will be inundated if a flat flood level is given. Lets say 1m flood from sea level, which is 0 m in this area. I was able to obtain/filter the polygons based on their elevation levels using query builders. However, this does not represent any sort of hydrological connectivity of these polygons. For an example, a low-lying area below 1m, surrounded by a ring dyke (dyke height above 1m, is also showing as inundated, with this method as this ignores hydrological connection and only consider elevation).
Is there any way that I can do this in QGIS and/or using others like GRASS or some tool? (except going for a dynamic hydrological/ hydrodynamic model like HecRAS)

Here, in the attached image, it is showing the polygons I have filtered from the whole domain, using elevation based filtering. Also, the rounded one is the problematic one, which indicates as inundated, but not hydrologically connected to other polygons. The middle white colour strip is a river stream, not included in the land area polygons