texts:
I am new to developement in QGIS. I use a the 3.4.4 version of QGIS.
My vector layer is like this:

So I want to know how to filter this data by date and apply the changes in the GUI. Here is what I tried but it is not working:
layers = QgsProject.instance().mapLayers()
for layer_id, layer in layers.items():
    idx= layer.fields().indexFromName("dateposte")
    values = layer.uniqueValues( idx )
    expr = QgsExpression("dateposte" LIKE '2019-02%')
    for f in layer.getFeatures():
        f[idx] = expr.evaluate( f )

EDIT:
I change this code and can find how many of the point pass the filter, the trouble is that there is no change in the map. So can anyone show me how to apply the filter in the vector on the map? Here is the new code
layers = QgsProject.instance().mapLayers()
for layer_id, layer in layers.items():
    if (layer.name() == "csv_points_acc") or (layer.name() == "csv_points_amb"):
        print(layer.name())
        expression = """"dateposte" like '%2019-01%'"""
        request = QgsFeatureRequest().setFilterExpression(expression)
        matches = 0
        for f in layer.getFeatures(request):
           matches += 1
        print(matches)


texts: ArcGIS: How to identify corner buildings

texts: In ArcGIS I have a polygone with buildings (from osm) and I have to select all the corner buildings. since I am new with GIS I don't know how to solve that. Tried to buffer all the intersect points of the streets, assuming that the the buildings within the buffer are all corner points, but to many non-corner-buildings are selected with that methode. Has anyone an idea how I can solve that?


texts: I am making a selection with ArcObjects, is there an Issue with iSelectionSet:AddList or do i have an error somewhere?

texts:
I collect a set of OIDs into an OID array (long), this happens in a loop and if i find an OID i want to add to the array, I REDIM PRESERVE to expand the array.

iOIDListCount = iOIDListCount + 1
ReDim Preserve iOIDList(iOIDListCount)
iOIDList(iOIDListCount - 1) = pRow.OID
Debug.Print(iOIDList(iOIDListCount - 1))

Once I use this array to make my selelction via iFeatureSelection.SelectionSet.AddList there is no error but a the resulting selection is completely wrong. Instead of getting an expected 89 records the attribute table is reporting 40 is selected but only about 10 records is physically selected on screen and in the attribute table (see figure 1)

NOT WORKING CORRECTLY:

pFeatureSelection.SelectionSet.AddList(iOIDListCount, iOIDList(0))

However when I use a For Next loop instead and add the OID's individually, I get exactly what I expect: 89 records. (see figure 2)

WORKS FINE:

    For i = 0 To iOIDListCount - 1
        pFeatureSelection.SelectionSet.Add(iOIDList(i))
    Next i

Figure 1 - Incorrect selection result when using AddList approach

Figure 2 - COrrect selection result when using Add (individually in a loop) approach


texts:
I have a bunch of 3D models drawn in Autocad MEP, which have the coordinates information but are not georeferenced.
I would like to be able to put them on a GIS map, so I can visualize both the lanscape and urban information I have, with the new Buildings I need to put on it.
Does anyone know how to drive the process?
It would be very helpful to find a replicable way since there are several 3Ds.


texts:
It generates the following error:

  línea 17: cvc-complex-type.2.4.a: Se ha encontrado contenido no válido
  a partir del elemento 'PropertyIsEqualTo'. Se esperaba uno de
  '{"http://www.opengis.net/ogc":spatialOps",
  http://www.opengis.net/ogc":comparisonOps ",
  http://www.opengis.net/ogc ": logicOps",
  http://www.opengis.net/ogc":FeatureId} '.

My code:

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0"
                       xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
                       xmlns="http://www.opengis.net/sld"
                       xmlns:ogc="http://www.opengis.net/ogc"
                       xmlns:xlink="http://www.w3.org/1999/xlink"
                       xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <NamedLayer>
    <Name>asia_cap_por</Name>
    <UserStyle>
      <Title>asia_cap_por</Title>
      <FeatureTypeStyle>
        <Rule>
          <Name>Baja</Name>
          <Title>Baja: 1.00 - 2.00 Kg/cm2</Title>
          <ogc:Filter>
            <PropertyIsEqualTo>
              <ogc:PropertyName>valor</ogc:PropertyName>
              <ogc:Literal>baja</ogc:Literal>
            </PropertyIsEqualTo>
          </ogc:Filter>
          <PolygonSymbolizer>
            <Fill>
              <CssParameter name="fill">#FFAA00</CssParameter>
            </Fill>
          </PolygonSymbolizer>
        </Rule>

        <Rule>
          <Name>Media: 2.00 - 3.00 Kg/cm2</Name>
          <Title>Media: 2.00 - 3.00 Kg/cm2</Title>
          <ogc:Filter>
              <PropertyIsEqualTo>
                <ogc:PropertyName>valor</ogc:PropertyName>
                <ogc:Literal>media</ogc:Literal>
              </PropertyIsEqualTo>

          </ogc:Filter>
          <PolygonSymbolizer>
            <Fill>
              <CssParameter name="fill">#FFFF73</CssParameter>
                         </Fill>
          </PolygonSymbolizer>
        </Rule>

        <Rule>
          <Name>Alta</Name>
          <Title>Alta: > 3.00 Kg/cm2</Title>
          <ogc:Filter>
            <PropertyIsEqualTo>
              <ogc:PropertyName>valor</ogc:PropertyName>
              <ogc:Literal>alta</ogc:Literal>
            </PropertyIsEqualTo>
          </ogc:Filter>
          <PolygonSymbolizer>
            <Fill>
              <CssParameter name="fill">#4CE600</CssParameter>
            </Fill>
          </PolygonSymbolizer>
        </Rule>

        <Rule>
          <TextSymbolizer>
         <Label>
           <ogc:PropertyName>descrip</ogc:PropertyName>
         </Label>
         <Font>
           <CssParameter name="font-family">Arial</CssParameter>
           <CssParameter name="font-size">10</CssParameter>
           <CssParameter name="font-style">normal</CssParameter>
           <CssParameter name="font-weight">bold</CssParameter>
         </Font>
         <Halo>
           <Radius>3</Radius>
           <Fill>
             <CssParameter name="fill">#FFFFFF</CssParameter>
           </Fill>
         </Halo>
       </TextSymbolizer>
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>


texts:
I've a list of a few hundred cities with their latitude/longitude.  Given another location (also in lat/long) I need to find the nearest city.

As I don't use any GIS, by now the obvious algorithm is to make a loop for all the cities, calculating the distance between the points.

Making the loop is practicable for me, but there is some easy to implement algorithm to accomplish that more efficiently? Or some light Java library that can help to solve that?

Notes: I don't need/want a complete GIS solution or heavy/complicated library. I prefer a less good but easiest and lighter solution because that is the only thing that I need to solve.


texts: Highlighting a specific country on a world map

texts:
I am using leaflet to display a list of country names in a side panel. I am also showing them highlighted on a world map.

When I click on a name in the side panel, I want to change the color of the associated country on the map.

Here is a snippet of my code:

DISPLAY LIST OF COUNTRY NAMES IN SIDE PANEL

$('#sidecol').empty();
$('#sidecol').html(COUNTRIES)
for (i=0; i < featureCollection.features.length; i++){
    //console.log(featureCollection.features[i].properties.NAME);
    $('#sidecol').append("div class = cntrys id = cntry_" + i + 'center>
' + featureCollection.features[i].properties.NAME + '/a>/center>br>')
cntryIds.push(featureCollection.features[i].id) //country id
}

COLOR FILL FOR COUNTRY CHOSEN

function getColor(d) {
cntryIds.forEach(function(value){
    //alert('index: ' + index);
    alert('D is ' + d);
    alert(value);
    return d == value ? 'green' :
                        'orange';
    });
}

STYLE FUNCTION CALLED WHEN COUNTRY NAME IS CLICKED ON IN SIDE PANEL

function newstyle(value) {
    //alert('1 ' + value);
    return {
        fillColor: getColor(value),
        weight: 2,
        opacity: 1,
        color: 'white',
        dashArray: '3',
        fillOpacity: 0.7
    };
}

CLICK EVENT OF COUNTRY NAME IN SIDE PANEL

$(".cntrys").each(function(index, value) {
  //console.log('cntrys' + index + ': ' + $(this).attr('id') + ' xxx:  ' + cntryIds[index]);
    $(this).click(function(){
        //alert($(this).attr('id') + " country id: " + cntryIds[index]);
        countries = L.geoJson(featureCollection,
        {
            style: newstyle(cntryIds[index]),
            onEachFeature: onEachFeature
        }).addTo(map);
        //console.log(getData());
        map.fitBounds(countries.getBounds().pad(0.5));
    });
});

I am trying to change the fillColor option in the newstyle function via the getColor function. Any assistance or suggestion is greatly associated...


texts: Why does my text symbol change when .mxd is converted to pdf?

texts:
I copied a small bit of text from a word document into a label which had the ' ← ' symbol in it (2190, alt+x).

Whenever I export the map to a pdf, the '←' symbol becomes a spiral. If I make the arrow symbol bold, it becomes a fleur de lis in the pdf. I have tried checking the 'Embed All Document Fonts', 'Convert Marker Symbols to Ploygons', changed all the Picture Symbol settings and nothing seems to keep the arrow from becoming something else.

It does not appear to matter if the text is a label taken from the attribute table, a Label Expression or is an annotation as the end result is the same.

If I save the word document as a .pdf, the arrow remains an arrow and does not change. Additionally, if the text with the spiral/fleur de lis is copied from the map pdf and pasted into Notepad, the spiral reverts to an arrow.

I realize there are many, many workarounds to this problem (not using the arrow, inserting a symbol etc), but I am more curious as to why it only seems to happen in ArcMap. It is now a matter of stubbornly trying to keep the arrow from being converted by ArcMap as a matter of principle.

Does anyone have any insight on why this occurs / how to keep the arrow?

Solution

By using

"<FNT name = 'Arial'>" +  ' ←' + "</FNT>" + ' Pathway to Moosehide Village '

it came out fine. Thanks for all the suggestions.


texts: Renaming columns in a QGIS attribute table with Python

texts:
I have some vector shape files with many columns.
I need to change the name of each column in the header from the old name to a more generic "Img1", "Img2", "Img3", etc progressively, starting from column #5 and leaving the name of columns 1 to 4 unchanged.

I wrote this Python code to be used with QGIS but it doesn't work.

layer = iface.activeLayer()
nImg = 1
nCol = len(layer.pendingFields())
cont = (5,nCol)
for field in cont:
    with edit (layer):
        oldName = layer.fieldNameIndex(field.name())
        newName = "Img" + str(nImg)
        layer.renameAttribute(oldName, newName)
    nImg = + 1

The error is

oldName = layer.fieldNameIndex(field.name())
AttributeError: 'int' object has no attribute 'name'

I tried also changing

with edit (layer)

to

layer.startEditing()

and unindenting the rows, but it doesn't work anyway.


texts: Spatial join between point files using FME

texts:
I have two files: the first has x,y,z information and is a (big) point cloud file, and second has only x,y information and is a (small) point shp file. The two files occupy the same area.

For each of the x,y points in the shp file, I would like to find the nearest x,y point in the point cloud file and pass the corresponding point cloud z value to the attribute table of the shp file. I would also like to avoid converting the point cloud to any other format (due to file sizes)

I am fairly new to FME and I'm not sure what is the right tool for the job. After some research, I tried the SpatialRelator but this seems to be the wrong tool (at least by itself) because it does not seem to support a "near" (within a certain distance) calculation. Maybe I'm too used to ArcGIS, but FME does not seem to have a "spatial join" type tool, or at least they are calling it something else...


texts: How to perform an anisotropic Elliptical Inverse Distance Weighting (EIDW) interpolation in R?

texts:
River bathymetry is generally more variable in the transverse than the streamwise direction, so I need to use anisotropic interpolation technique like Elliptical Inverse Distance Weighting (EIDW) (see Anisotropic considerations while interpolating river channel bathymetry) to construct a digital elevation model (DEM) of the surveyed channel. I only found how to perform an Inverse Distance Weighted using the gstat package in R:

library(sp)
data(meuse)
coordinates(meuse) = ~x+y
data(meuse.grid)
gridded(meuse.grid) = ~x+y

# inverse distance weighted:
x <- krige(log(zinc)~1, meuse, meuse.grid)


texts: Set up Esri Proxy on Visual Studio for Web Appbuilder

texts: I am trying to set up Esri proxy for my Web Appbuilder app on Visual Studio , anyone has done that? Setting up proxy tutorials show setting up on IIS7, VS runs IIS Express, therefore the set up would be different


texts:
I'm facing a strange behavior using buffers in Turf.js.

When trying to do a 50m buffer around a point, the result is actually closer to 35m.

On the screenshot below, the smaller circle is a turf.buffer() and the bigger one is a turf.circle() using the same radius/distance. This second one is correct.

Here's my code :

var map = new ol.Map({
   target: 'map',
   layers: [
     new ol.layer.Tile({
       source: new ol.source.OSM()
     })
   ],
   view: new ol.View({
     center: ol.proj.transform([5, 45], 'EPSG:4326', 'EPSG:3857'),
     zoom: 17
   })
});

var distance = 50

var center = turf.point([5, 45])

var bufferedWgs84 = turf.buffer(center, distance/1000, {units: 'kilometers'});

var circleWgs84 = turf.circle(center, distance/1000, {units: 'kilometers'});

var buffer = turf.toMercator(bufferedWgs84)

var circle = turf.toMercator(circleWgs84)

var layer = new ol.layer.Vector({
                source: new ol.source.Vector(),
        style: new ol.style.Style({
                fill: new ol.style.Fill({
                color: 'rgba(255, 255, 255, 0)'
                }),
                stroke: new ol.style.Stroke({
                color: '#737373',
                width: 2
                })
            })
            });
map.addLayer(layer);
layer.getSource().addFeature(new ol.format.GeoJSON().readFeature(buffer))
layer.getSource().addFeature(new ol.format.GeoJSON().readFeature(circle))

Fiddle : https://jsfiddle.net/fhkgzd2o/1/


texts: Merging features into single programmatically

texts:
I have a shapefile with polygons.
I can merge them manually by going into edit session and, selecting a some features to be merged, and then using merge button under Editor.

How to do this programmatically (in VBA)?

I've found the analogous ways but it doesn't solve the problem:

Standard Dissolve tool (from ArcToolbox) doesn't work correctly for this purpose, because it cuts off many fields and because it operates with whole FeatureClass but not with a separate features.
As I know by reading some examples, using an IBasicGeoprocessor (VBA)  have the same problem - it operates with whole FeatureClass, producing a new FeatureClass, but not with separate features, mearging them within the same FeatureClass.

@kenbuja, Here is my try in vba (not in vb.net) :
   Using a CType causes the Compile error: Sub or Function not defined.
   What do I do wrong?

Dim pMxDoc As IMxDocument
Dim pCountyLayer As IFeatureSelection
Dim pCountySelection As ISelectionSet
Dim pCountyCursor As IFeatureCursor
Dim pCountyCursor2 As IFeatureCursor
Dim pFirstFeature As IFeature
Dim pSecondFeature As IFeature
Dim pThirdFeature As IFeature

Set pMxDoc = ThisDocument
Set pCountyLayer = pMxDoc.FocusMap.Layer(1)
Set pCountySelection = pCountyLayer.SelectionSet

' Create a cursors from the selected feature in Layer (1).
pCountySelection.Search Nothing, True, pCountyCursor
pCountySelection.Search Nothing, True, pCountyCursor2
' Return the selected features from the cursor.

Set pFirstFeature = pCountyCursor.NextFeature

Set pThirdFeature = pCountyCursor2.NextFeature
Set pSecondFeature = pCountyCursor2.NextFeature

Dim pFirstGeometry As IGeometry
Set pFirstGeometry = CType(pFirstFeature.Shape, IGeometry)
Dim pTopoOperator As ITopologicalOperator
Set pTopoOperator = CType(pFirstGeometry, ITopologicalOperator)
Dim pSecondGeom As IGeometry
Set pSecondGeom = CType(pSecondFeature.Shape, IGeometry)
Dim pMergedGeom As IGeometry
Set pMergedGeom = pTopoOperator.Union(pSecondGeom)


texts: Delete rows in which the cell of a specific column is empty

texts:
Is there a way to delete all rows in QGIS in which the cell of a specific column in the attribute table is empty? e.g. line three in the table below. I need to do this for a large data set with about 180,000 rows.

| ID  | 2LE  | CHECK |
|-----|------|-------|
| 1   | PT   | FALSE |
| 2   | PT   | TRUE  |
| 3   | NULL | TRUE  |
| 4   | AT   | FALSE |


texts: Field added to SpatiaLite table doesn't show in QGIS

texts:
I have a spatial table in a SpatiaLite database that I'm editing in QGIS with no problem (the spatial table was created through SpatiaLite GUI).

I then added a text field to the table through SQL statement and I can edit this column in SQLite Studio or SpatiaLite GUI! However, this column doesn't show in QGIS, even if I add the layer to a clean project!

Any thoughts on how to update the field list in order to be correctly recognized in QGIS?


texts: Sorting layers in legend using python in QGIS 2.4

texts:
In QGIS 2.2 , I was able to sort in alphabetical order the layers simply using this:

self.maplegendtree = iface.mainWindow().findChild(QTreeWidget, "theMapLegend")
self.maplegendtree.sortItems(0, Qt.AscendingOrder)

I search a lot and I understand there where some changes with the legend in QGIS 2.4 ...

I tried this :

self.maplegendtree = self.iface.mainWindow().findChild( QTreeView, 'theLayerTreeView')
self.maplegendtree.setSortingEnabled(True);
self.maplegendtree.sortByColumn(0, Qt.AscendingOrder);

No error ... but no changes in the legend or mapcanvas?


texts: How to convert ArcMap color ramps to QGIS?

texts:
I have both ArcGIS 10.0 and QGIS 2.0

How can I reproduce a specific ArcGIS color ramp in QGIS?

If possible, I would like to simply convert the Esri *.style file but perhaps there is another approach that I am not aware of.

I would prefer not to have to manually create a new QGIS color ramp from scratch.


texts: ArcGIS Server 10.2.2 - Map Export Sometimes Returns Blank Image

texts:
I've run into a fun little issue with ArcGIS Server 10.2.2 to end my week with...  I have discovered today that the export operation (the one that renders map layer images) is randomly returning blank images.  As in - I enter a export url into the browser, hit F5, get expected image.  Hit F5 again - get blank image.  Hit F5 again - get expected image.  None of the parameters are changing - its exact same url.  There is no observable pattern - sometimes I got to hit F5 multiple times before I get expected image.

I remember facing a similar bug back at 10.1 - where every other GET query failed (How does that make it through testing???).  This seems awfully similar to that.  I starting search GeoNet - and have seen similar posts.  Does anyone know if there is any official NIM for this issue - and more importantly a solution or workaround?


texts: How do GPS/Glonass/Galileo satellite know their position?

texts: We use satellites to calculate our position. How do satellites know their position if it is right or wrong ? I have heard that the base station records and corrects the satellite path but how do they know what is right path what is wrong path ?


texts:
I get this error when undertaking in query in QGIS2.4 or QGIS2.6:

An error occurred when executing the query.
The data provider said:
OGR[3] error 1: Unrecognised field name Fire.isochrones_2002_BigDesert_21122002.OBJECTID.

All my field names are less than or equal to 8 characters.

Anyone got any ideas on how to fix this?


texts: Database has duplicate values after running pgrouting

texts:
I am a newbie with postgresql and sql statements so maybe it is a trivial question, I'm sorry for that.

I am working on a postgis database and try to use the pgrouting function.  My goal is to store the results of the shortest path query in a different table than the original one, which is working so far. The problem is, that it keeps adding the results to the original table too, so there are duplicate values in it afterwards, which I try to avoid but without success so far. I don't know if it should be like this? Anyway can you suggest any solutions to this problem? (I solved it so far so that I delete the duplicate values afterwards, but I really want my database to stay unedited... is it even possible?)

I am using postgresql 8.4 with postgis 1.5 and pgrouting 1.05

This is my query:

 INSERT INTO
route (

SELECT
topo.id,
topo.osm_id,
topo.osm_name,
topo.osm_source_id,
topo.osm_target_id,
topo.clazz,
topo.flags,
topo.source,
topo.target,
topo.km,
topo.kmh,
topo.cost,
topo.reverse_cost,
topo.x1,
topo.y1,
topo.x2,
topo.y2,
topo.geom_way,
topo.available

FROM
Shortest_path('select
*
from
testro_2po_4pgr',
2345,
9234,
false,
true)
AS path,
testro_2po_4pgr
AS
topo
WHERE
path.edge_id=topo.id);

Thank you very much in advance!


texts: PostGIS: how to split a polyline into straight segments?

texts:
I'm working with PostGIS 2.1.8. I need to split a polyline into as many straight segments this polyline contains. Here is an illustration to make it clear enough:

The initial feature is only one row and the query should produce a 4 rows result with attributes values being maintained and geometry column changed (I don't want to update length column):

# Initial table:

id   | name    | length   | the_geom
-----+---------+----------+-----------------
1    | Line 1  | 12.2     | 01000111100000

# Result table:

id   | name    | length   | the_geom
-----+---------+----------+-----------------
1    | Line 1  | 12.2     | 01000111100000
2    | Line 1  | 12.2     | 01001111000011
3    | Line 1  | 12.2     | 01000110010000
4    | Line 1  | 12.2     | 00000111100000

I have found the function st_subdivide() (see here) but it is available from PostGIS version 2.2... I can't upgrade my PostGIS version because I'm at work and not allowed to do it.

Any idea on how to achieve that with an older than 2.2 PostGIS version?


texts: Count overlapping polygons in single shapefile with R

texts:
I have a shapefile of species ranges for ~ 346 species.  Some species have multiple polygon ranges so there are 1,078 polygons and many polygons are multi-part.  It looks something like this:

I would like to get the max count of overlapping polygons anywhere within the layer bounds. There are two ways I have considered accomplishing this, but I wanted to see if there is a simple way before I undertake the effort to try either.  I would prefer to work in R, but could potentially use ArcGIS 10.2.

I tried using the 'Count Overlapping Polygons' toolbox from http://www.arcgis.com/home/item.html?id=1dd4a6832b3d40b494dbf8521cc5134c which seems like it would allow me to extract the information I need.  However, it failed with an 'Out of memory' error.  I have a Winx64 computer with 8GB of RAM and only had a few other processes running.

I came across this post by @Roger Bivand but it appears to only function between two shapefiles, which will not work for me.

The other approach that I thought might work is to rasterize each polygon and then sum all the rasters.  I found a function here that allows the conversion to raster.  I thought this might be a very time-consuming and computationally expensive, so I have not yet tried that approach.


texts: How to change length to 2 decimal places?

texts: I'm trying to work out how to calculate the length of a line to only 2 decimal places in QGIS. When I use the field calculator to create a new field, it doesn't matter what changes I make to the field type, field width or precision, I continually get the length to 13 decimal place.


texts: Split up shapefile into little shapefiles

texts:
I have a very large shapefile (all the roads in RI). I need to split it up into chunks so that I can process the file with multiple worker threads in Python. So far I have created a grid in QGIS that is overlaid on the vector data. I then ran the intersection algorithm with the grid and the shapefile so that the vectors are split up and are assigned an 'id' for the box they're part of.

Now, since all the vectors in a particular box of the grid all have the same 'id', can I extract them and save them to their own shapefile automatically? I am aware of selecting all vectors with an 'id' and then copying them to their own layer and then saving that layer, but I have so many boxes/id's that it would take forever to save each of them off. How can I use QGIS to automatically make a new shapefile for each id/box?

Basically I want create a new shapefile for each distinct value of a particular field in the shapefile attributes.

Original Shapefile

Shapefile after grid intersection

A set of vectors which share the same bbox/id


texts:
I'm trying to perform a simple buffer on GeoJSON file using Python script. I am able to do it on a shapefile using Python script and on GeoJSON using QGIS fix distance buffer tool but I can't figure out how to do it on GeoJSON using Python script. I read many posts on it but none of them include working solely on GeoJSON. I also know that it is possible to convert GeoJSON to .shp but I need to work only with GeoJSON files.

My code that works for .shp file:

from osgeo import ogr

ds = ogr.Open(r'C:\vector files\tryyy.shp', 1)
in_lyr = ds.GetLayer()
out_lyr = ds.CreateLayer('new_buff44',in_lyr.GetSpatialRef(),ogr.wkbPolygon)
out_lyr.CreateFields(in_lyr.schema)
out_defn = out_lyr.GetLayerDefn()
out_feat = ogr.Feature(out_defn)
bufferDist = -0.00015075

for in_feat in in_lyr:

    geom = in_feat.geometry()
    out_feat.SetGeometry(geom)
    for in_feat in in_lyr:

        geom = in_feat.geometry().Buffer(bufferDist)
        out_feat.SetGeometry(geom)
        for i in range(in_feat.GetFieldCount()):
            value = in_feat.GetField(i)
            out_feat.SetField(i, value)
        out_lyr.CreateFeature(out_feat)

del ds
print "finish"

When I change the file from shp to geojson in the org.Open line I get:

  NoneType' object has no attribute 'CreateFields'

I thought that if I read vector file as OGR datasource object then it doesn't matter which format is the original vector file.

Can someone help me understand what I should change so that the code will work on GeoJSON files?


texts: Change cluster label on Leaflet web map generated using QGIS2WEB plugin

texts:
I want to change cluster labels on Leaflet maps generated using QGIS2WEB plugin. Actually I want these labels to show a sum of values from point attribute fields (all points on same location have that field, but different values). Also, if that is possible I want to change cluster colors depending on sum value (bigger values shades of red or orange colour, lower values shades of yellow or green) some kind of heat map.
My question is what I should change in generated code to solve my problem?


texts: Model to query all feature classes in geodatabase by selecting by attribute?

texts: I have a geodatabase that has 60+ feature classes that each contain records for 60+ locations.  I need to export all records from each feature class for a specific location.  In each FC there is a common attribute that identifies each location (i.e. instln_id).  I would like to do this in ModelBuilder but have not had any luck with iterating feature classes or workspaces and selecting by attributes.  Basically I want to scan through every feature class for records with a specific attribute and export records to new feature classes.  Does anyone have any ideas or examples of something you've done that is similar?


texts:
I hope you can help me with the following:

We need to create a stand alone asp.net based site that will display some maps. We are using IIS Express to run a site and now we need to add the map server part.

We have map files from the existing Apache/php & OpenLayers based site and need to add them to our site.

My MapServer understanding is basic so any help is appreciated.

Also where I can find any good references/tutorials?

Thanks in advance.


texts: Raster saving symbology break values manually in ArcGIS Desktop?

texts:
I am trying to manualy change the break values for my raster image through layer properties and then symbology and classified manual. Whilst i am able to do this, I cannot figure out a way to save these values so that the next time I am on ArcMap it loads the raster with the manually changed break values.

Does anybody know how to do this?


texts: Creating map gores in ArcGIS for Desktop?

texts:
Does anyone have clear instructions on creating gores for a globe in ArcGIS?

I have tried to work with these:

http://mappingcenter.esri.com/index.cfm?fa=ask.answers&q=1725

But have become lost. I have been using this world border shapefile:

http://thematicmapping.org/downloads/world_borders.php

Can anyone give some clear advice on how to create a set of globe gores in ArcGIS?


texts: Using ArcPy and certificate-based login?

texts:
We have an ArcGIS Server/Portal federation, and I'm trying to publish a map using Server.  The login gauntlet is that Server defers to Portal, which, in turn, defers to a SAML server which checks the certificate presented and makes the go/no-go call.

The certificate needed was imported by the user.

When I run my publish script, it runs a while, then pops up a dialog box that basically says, 'Is this the proper identity certificate?' with enough info to make that determination, with OK/Cancel options.

That's fine for current testing, but not so great for when we want this process to be run in a fully-automated fashion.

So the question is, is there some way, in arcpy, to basically predeclare, 'This file is the certificate you want, present it if challenged'?


texts: Creating reclass-rules-textfile within Simple Python Editor of GRASS GIS

texts:
We are writing a script with the simple Python editor (GRASS GIS).
We try to set the rules for a reclass command within the script (without creating a txt.file manually).

Is it possible to set the rules within the script or is it possible to create a txt.file with the simple Python editor?


texts: Calculating distance between set of points: lat and long in R?

texts:
I have shapefiles containing the location of several individuals along several months. Positions were documented with a 30 minutes interval.
 Tables contain the following data fields: collar (ID), date, time, lat, long, altitude.

I would like to compute the distance between each pair of points of each individual, which have been taken by 30 min or 1 hour interval apart?
a number of positions and want to calculate the length between the positions - per hour, and per day.

Can someone guide the way to achieve it using R?


texts: Selecting first items of same attribute value using PyQGIS?

texts:
I've been getting into PyQGIS, and as a start I want to select the first instance of pair values of a point layer (start and end of a street), identified by the same integer, I got this, but it doesn't work:

("int_id" is the value that is the same for the start and end points)

cap = iface.activeLayer()
feat = cap.getFeatures()

frm = []

for f in feat:
    if f["int_id"] not in frm:
        frm.append(f.id())

cap.select(frm)

I need to create lines from these points (start and end), I haven't found a way to create them in PyQGIS, but I found point connector plugin that could at least create the streets that are a straight line, but first I need to separate the start points from the end points.


texts:
When I connect to my postgis database and edit/add features and keep the edit mode open, many times it loses the connection.

I cannot save my new/edited features and I need to restart the software and start from the beginning. I tried many type of qgis versions and only the 2.12 qgis lyon able to reconnect by clicking on the elephant head and connect and after that I can save the database. With the other versions I just save after every feature and hope it will save the next one as well. Is there any qgis settings or database settings which allows the connection to keep open? We are using qgis on multiple computers and all of them have the same problem... this question is relating to this one but no solution:  QGIS doesn't handle PostGis connection resets?

this is the error generated:


texts: Using ArcMap to find borehole location to drill bore to fetch underground water?

texts: I have got an raster img file. I have the corresponding KML file. I wanted to find underground water using these two files. I wanted to know the steps to identify the location of water.


texts: Finding multiple least-cost-paths for a network of several points using QGIS?

texts:
How can I find the least-cost-paths between a network of several points using QGIS or other easily available tools? That is, the easiest way to get from each point to every other point in the network, as opposed to just between one start and one end point.

The start and end points are polygons scattered unevenly across the map, and the intervening landscape is a patchwork of fields etc. for which the permeability score of each landscape feature is known - the least-cost-paths will need to choose any possible route through the fields, not follow roads or similar.

(I am using QGIS-2.14.1-Essen with GRASS).


texts: Classifying DEM base on lowest elevation of stream generated by flow accumulation of same DEM using ArcMap

texts:
How do I classify a DEM base on the lowest elevation of a stream (stream generated by flow accumulation of same DEM)?
I want to find those elevation of DEM that are close to the floor of the stream. It's not a simple classification because the stream floor does not have same elevations in its whole lenght, meaning that there is a big elevation difference between start and end point of stream (Stream has alot of slope!)
Stream can be a river or a wash.


texts:
According to the API:

ol.extent.containsExtent(extent1, extent2)

checks to see whether extent1 is contained by extent2.

With the following variables:

extent1 = [-93, 44, -94, 45];

extent2 = [-95, 43, -92, 46];

ol.extent.containsExtent(extent1, extent2) returns false, though extent1 is clearly within extent2. What am I doing wrong?


texts: Rasterize will output 0 as NODATA (QGIS)

texts: I have a data range which includes 0 ex +3 to -4. I would like the rasterize to include my 0 data (create a pixel) and just exclude the NODATA (where there is no points)


texts: I'm trying to get QGIS (1.8) to scale the point markers for each town in a district proportionately to the size of their population (i.e. bigger population = bigger circle). I can get them to vary by colour by graduating the properties for the layer with the population data. When I try to use Advanced>size scale field function though I either get nothing or the entire screen goes one colour. I've not found a way in QGIS to scale by size only which is what I really want to do. Ideally I want the circles to graduate in size only and to have the same colour. I'm hoping that there is a straightforward way of doing this in QGIS.


texts: Opening ArcGIS Desktop on specific license level?

texts:
We are on Citrix and would like to provide 3 shortcuts to our users so they can launch the three license levels (ArcView, ArcEditor and ArcInfo).  We can't do this through desktop admin as the users don't have admin rights and we would like an intuitive workflow anyway.

How do others provide access to the different license levels in their organisations for the users to quickly open the product license that they want?


texts: Action on group of features from layer in QGIS?

texts:
I need to implement a rename action on a selected group of features in a certain order. The idea is that I would rename one feature and then select the features that are 'next in line' and execute a rename function (if i name the first one A then the next would be B, C, D, etc). My problem is that a python layer action is run on a single feature and I'm not sure that the feature sequence in which the action is done is also the proper one. I would like to avoid as much user input as possible.

Could you advise me the best way to do this? I need to distribute the solution to my user group.

QGIS version 2.8.1


texts:
LAStools program only limited to UTM platform for coordinate transformation
(las2project, las2transform).

Is there any other ways to transform coordinate for local projection in LAStools? For example RSO Malaysia.


texts:
I am trying to extract current USGS DEM data that corresponds to an old USGS map. I have downloaded the DEM covering the corresponding area, and am now trying to extract the correct subset. The area covered in the old map is irregularly shaped. I have successfully extracted a rectangular area using GDAL, however this includes extra data not present in the original map.
My understanding is that GDAL translate only accepts rectangles as bounding boxes.
Is there a way to set all the values in the unneeded region either to a no-data value or to 0?
Here is a picture illustrating the issue:


texts:
I am looking to replicate a tool that is available in Geospatial Modelling Environment (GME) within QGIS, preferably not using GRASS. The tool is isectlinerst (Intersect Lines With Raster) http://www.spatialecology.com/gme/isectlinerst.htm

I am hoping to have a tool that can extract raster information (i.e. slope) for an entire linear network (i.e. roads). I have seen other post on here that show a rasterize process, but I am looking to maintain the vector input data and add the attributes to the data.


texts: Find max and min values from a table based on two fields

texts:
I have a geodatabase table with the following fields: County_Name, State_Name, Date_Time, Customers

Example

Alameda, CA, 8/1/2014 12AM, 6
Alameda, CA, 8/1/2014 1AM, 1
Alameda, CA, 8/1/2014 2AM, 3
....

Data for all US counties is in the table. I just want to create a new table that contains the

county, state, max customers for the day (24 hour period), time of the max customers
county, state, min customers for the day (24 hour period), time of the min customers

Many states have the same county names as other states.

Do I need to create a cursor to obtain a list of the unique county and state names first?

This will give me a list of the county names but I also need the state names.

with arcpy.da.SearchCursor(tablesMerged, ["county_nam"]) as cursor:
        print sorted({row[0] for rows in cursor})


texts:
I am developing a plugin for QGIS 3.10.
In part of my plugin, I have a UI using PyQt.
This form needs to be a singleton.
But when I have click several times on the plugin it generates several forms
Here is my code for UI
class AddDuctUI(QtWidgets.QDialog, FORM_CLASS):
    __instance = None
    def __new__(cls,*val):
        if AddDuctUI.__instance is None:
            AddDuctUI.__instance = QtWidgets.QDialog.__new__(cls)
        AddDuctUI.__instance.val = val
        return AddDuctUI.__instance

    def __init__(self, parent=None):
        """Constructor."""
        super(AddDuctUI, self).__init__(parent)
        self.setupUi(self)

When I check the object it is single instance but I have several windows shown.
This is the way I call the UI:
def show_form(self):
     main_dlg=AddDuctUI()
        self.main_dlg.show()
        result = self.main_dlg.exec_()

        if not result:
            self.reset()

        if result:
           self.apply_changes()
           self.reset()

How can I make sure only one window from this UI is shown at the time?


texts: Info Window close button color change in CartoDB?

texts:
I was looking at changing the infoWindow close button color -

a href="#close" class="cartodb-popup-close-button close" style="background-color: #000; font-color: #000;">x

The 'background-color' works but every attempt at changing the color of the 'x' has failed.
I was thinking the x might be a font but doesn't seem to be. font-color doesn't change anything.

Any thoughts?


texts: Is there a simple command line or python based way to visually register a pair of overlapping stereo images and clip the non-overlapping areas out?  These images do not have spatial information associated.  Thanks!


texts: How to merge lines from shapefile using PyShp and Shapely and save them in a new shapefile?

texts:
I have files with line segments that cover road network of towns. The issue is that each street with specific name consists of smaller segments that I need to merge into one line covering the whole street.
What I tried in the code below is to read a shapefile, create dictionaries where under a streetId are all of the street's segments and then merge it using ops.linemerge. The problem I have is when saving to the file in the following line
w.line(merged_line)

The error shown is
    w.line(merged_line)
  File "C:\OSGeo4W64\apps\Python37\lib\site-packages\shapefile.py", line 1977, in line
    self._shapeparts(parts=lines, shapeType=shapeType)
  File "C:\OSGeo4W64\apps\Python37\lib\site-packages\shapefile.py", line 2062, in _shapeparts
    for part in parts:
TypeError: 'LineString' object is not iterable

What object should I provide to the writer's line method? Is there any better way to merge the segments into streets?
import shapefile as shp
from shapely import geometry, ops
from shapely.geometry import shape

sf = shp.Reader("shapefiles/town1")
streetShapeDictionary = {}
streetSegmentDictionary = {}
shapes = sf.shapes()
records = sf.records()
fields = sf.fields

recordIndex = 0;
for record in records:
    streetId = record[0]
feature = sf.shapeRecords()[recordIndex].shape.__geo_interface__
linestringShape = shape(feature)
if streetId in streetSegmentDictionary:#add street segment to the list
    streetSegmentDictionary[streetId].append(record)
    streetShapeDictionary[streetId].append(linestringShape)
else:#add new key to the street dictionary
    streetSegmentDictionary[streetId] = [record]
    streetShapeDictionary[streetId] = [linestringShape]
recordIndex = recordIndex + 1

w = shp.Writer('shapefiles/whole_streets')
w.field('name', 'C')
streetShapesArray = []
for street in streetSegmentDictionary:
    for shape in streetShapeDictionary[street]:
        streetShapesArray.append(shape)
    multiline = geometry.MultiLineString([*streetShapesArray])
    merged_line = ops.linemerge(multiline)
    w.line(merged_line)
    w.record('streetName1')
    streetShapesArray = []

w.close()


texts:
I'm looking to merge two road layers. I am interested in the traffic volume in the second layer. The roads do not lay exactly on top of each other.

The issue is that sometimes when there is an intersection the road is taking the perpendicular value but I want it to take the parallel value. My first thought was to use a spatial join and join the closest join feature to the target feature with a 20 foot tolerance. I have also tried using the buffer tool then joining roads that are contained. But both methods are not producing the desired results.

An example photo is below. The buffered polyline is taking the value of the north and south road instead of the east and west road. I've also looked into Merge Divided Roads and snapping line to line but don't have the license for them.

I have ArcGIS 10.5 and a basic license


texts: ogr2ogr: encoding issues when importing from Oracle Spatial to PostGIS

texts:
I'm trying to import some spatial tables from Oracle Spatial to PostGIS using ogr2ogr, but I can't figure out why the special characters on the attribute tables get messed up while I'm doing it.

My environment is:

Oracle Database 10g / Oracle Spatial 10.2.0.4 NLS_CHARACTERSET WE8ISO8859P1, NLS_LANGUAGE AMERICAN
PostgreSQL  9.5 / PostGIS 2.2 installed on Windows 8.1 Pro PT-BR. Database encoding is UTF8, character type IS Portuguese_Brazil.1252

My ogr2ogr command is:

ogr2ogr -progress -f "PostgreSQL" -nln schema.table -skipfailures -overwrite -lco FID=ID_COLUMN -lco GEOMETRY_NAME=geom PG:"host=localhost user=postgres password=yadayada dbname=database" OCI:username/password@server:port/database:schema.table

The table is imported perfectly, except for the special characters like "~". All the "ã" are replaced by "?", for example. I've already tried to use SET PGCLIENT ENCODING=LATIN 1 and SET PGCLIENTENCODING=WIN1252 without success.

I don't know if it helps, but if I try to open the tables directly from Oracle into QGIS the characters are displayed properly.


texts:
I have run a linear regression and obtained an equation that describes the link between a group size and the area used by this group.
Now, I would like to enter this equation into the GIS to obtain buffers of different sizes depending on the group size. For instance, a group of 50 individuals may end up with a buffer of 100m diameter. I was going to add a field in the attribute table "buffer" and enter the equation.
After that, I was wondering if I could use my different buffers to extract values from rasters in a simple way?
I mean that I would like to extract the average value but based on my buffer size and without having to make a hundred shapes...
I hope that I have been clear

Thank you so much for your time,


texts: Automating map export from list of inputs using ArcMap

texts:
Does anybody know of a good premade ArcMap tool to iterate over a list of inputs (rasters, shapefiles, etc) and export to PDF/PNG?
Preferably a Toolbox tool that uses the current mxd as a template and possibly allows for symbology and/or legend options.
I have a routine that produces 50 raster surfaces, I just want to be able to flip through them without exporting 1 by 1.
The following had some good leads but seemed out of date:

Automating map creation for 50,000 unique routes using ArcGIS Desktop?
Automating entire map production in ArcMap using ArcPy/ArcObjects?


texts: Store human readable layer name in QGIS style files?

texts:
QGIS 2 and 3 as well seem not to store the Layer name (set via Layer Properties | Source | Settings | Layer name or via 'Rename' in the layers context menu) in styles files. I.e. if i have a PostGIS table 'flst_teil' and would like to have it in the Layers Panel in a nice human-readable form like 'Flurstück Teil'.

For QGIS 3 that astounds me even more, since herein I can explicitly select which configuration tab of a layer I wish to store:

Is this intented behaviour or a bug? Do I have a chance to achieve this or do I have to rename the layers any time I load them?


texts: Check if coordinates intersect rotated rectangle

texts:
In a non-GIS system, I have a Python script that checks to see if coordinates are valid:
def is_latlong_valid(laty, longx):
    #Verify if the numbers are legitimate UTM coordinates
    return (4856762 <= laty <= 4879371 and
            663546 <= longx <= 677126)

The script checks if the coordinates are within a certain un-rotated rectangle (grey text, black box).

I suppose this could be considered a sort of minimum bounding envelope.

In reality, the geographic area I want to check the coordinates against is  a rotated rectangle (black text, blue box).
Is there a way to enhance the script so that it checks the rotated rectangle, instead of just the un-rotated rectangle?

I do not have access to any spatial libraries.
I don't even have access to most of the standard Python libraries, and I am unable to import additional libraries.

Related questions:
Stack Overflow: Determine if points are within a rotated rectangle (standard Python 2.7 library only)
Code Review: Determine if points are within a rotated rectangle


texts: SAGA not working under QGIS 2.18.10 (OSGeo4W)

texts:
Recently, saga algorithm disappeared from the processing panel. However, they are still listed under processing->options.

How to make them available again?


texts: Standard colour ramp for wave height visualisation

texts:
I'm working on a visualisation of weather, including wave heights. The source data is from a WaveWatch 3 model, and I'm using GeoServer 2.7 (with NetCDF and dynamic colormap extensions) to do the rendering. My question relates to how to best represent that data to users who have domain expertise (e.g. commercial shipping / deep water fishing / Navy) but not much numerical forecasting expertise.

The use of the data isn't (completely) known, and will vary by user. One example is that training for helicopters or aircraft might not be allowed above seas that are above some (specified) height. Alternatively, it might be some relative indication with a visualisation over time (e.g. so fishing can occur in the smoother water to the west, before moving to the east after that area settles).

The end result will be WMS, which will be rendered above the basemap (with opacity that is user selectable from 0 to 100%). The basemap is user selectable on the client side and I have no way to know what the basemap colour will be (best guess: some solid colour for deep ocean, perhaps blue, white or black).

There is a good question (and some excellent answers) on an appropriate colour-ramp for DEM, which was informative. However most of the answers for that were related to making the terrain "look right" (e.g. with hill shading), which isn't really appropriate for something like wave height.

My priority is to communicate the weather conditions first, and to look good (as a close) second.

My question is: what is an appropriate wave height colour ramp? If there was an accepted international standard (or some fairly common national standard) I would likely prefer that. Failing that, I'd go with accepted designs for this kind of visualisation.


texts:
I am trying to place points in QGIS that are given in a csv-file, coordinates are given in WGS84 (4326) but the positions are like this
LAT;LONG
5724.3180;1548.4021

I can't get QGIS to place them correctly. When I divide them by 100 I get the extent correct but the place is not correct 57.24318 / 15.4021
It should be Hultsfred. The points represent trees along Väg 23 between Målilla-Hultsfred
How do I solve this?


texts: I am use AcGIS 10.1 and have a DTM with an 8 m grid resolution. Which tool could I apply to smooth the surface?


texts: How do I maintain the current map extent of a project file on save in QGIS2?

texts:
I'm using QGIS 2.0.1 on a Mac and having an issue with maintaining the map extent after saving the project.  I save my project file with a certain extent and when I re-open the file the map is centered on 0,0.  I can use the "Zoom Last" button to get back to the extent of the map when I saved.  This occurs regardless of how Enable 'on the fly' CRS transformation is set.  Layer order, symbol changes, etc are all maintained on save, only the extent jumps to 0,0.  My test shapefile is in WGS84.

Looking at the project file XML I see the correct extent coordinates in the <extent> tags.  So it looks like QGIS knows the correct saved extent (in the project file and "Zoom Last"), but jumps to 0,0 on re-load.

Thoughts on what could be going on here?  Thanks!

Update:
I opened the same project file on another machine (Windows) and the extent behaved as you would expect.  My Mac install of QGIS seems to be the issue here.  It is on OSX Mavericks and so that might be at the root of this.

Update 2:
On opening a new dataset with QGIS (right click on shapefile -> open with QGIS), the map jumps to center on 0,0 with an option to "Zoom Last" back to the extent of the new data.  Behaves as above, recognizes the initial extent and then immediately zooms to 0,0.


texts: How to insert a GeoJSON polygon into a PostGIS table?

texts:
I need to insert a polygon from GeoJSON to my PostGIS table. This is how the SQL query looks like.

INSERT INTO tablename (name, polygon)
VALUES (
    'Name',
    ST_GeomFromGeoJSON(
        '{
            "type": "Polygon",
            "coordinates": [
                [7.734375,51.835777520452],
                [3.8671875,48.341646172375],
                [7.20703125,43.580390855608],
                [18.6328125,43.834526782237],
                [17.9296875,50.289339253292],
                [13.7109375,54.059387886624],
                [7.734375,51.835777520452]
            ]
        }'
    )
)

Unfortunately, I get an error message.

ERROR:  Geometry SRID (0) does not match column SRID (3857)

The GeoJSON is already in the right reference system. But this isn't specified. How do I specify the SRID in the GeoJSON? What does the GeoJSON need to look like?

Update: When I wrap the geometry created by ST_GeomFromGeoJSON with ST_SetSRID(..., 3857) it throws another error. In my view it doesn't seem that the geometry has a Z dimension.

ERROR:  Geometry has Z dimension but column does not


texts: Removing datapoints from LAS dataset using ArcMap?

texts:
I'm working with lasd data sets in ArcMap 10.5 and I'm trying to create an density raster from the ground points.

How would I either remove all but the points classed as ground from the las data set or export the datapoints as a vector point layer?


texts: QGIS 3.4 Vertex tool highlight colour change

texts:
I can change the colour for selected items fine. I have even figured out how to change the colour of the preview when moving a vertex. How can I change the colour of the highlight when hovering over a feature with the vertex tool?
The default is red and the line features I most commonly work with are also red so I was really hoping there was a way to change this to make seeing the vertex I'm about to select easier.


texts: Arcpy RasterToNumPyArray does not read floating point raster

texts:
I have below rasters:

One with dimension (13646, 7510) where INTEGER pixel values ranges from 0-14.
One with dimension (13646, 7510) where DOUBLE (upto 10 decimal places) pixel values ranges from 0-181257.

I can easily read 1 with RasterToNumPyArray but when I try to read second raster exception raised as below:

  return _RasterToNumPyArray(*args, **kwargs)

  MemoryError:

I tried with scipy.misc it also gives me error for 2 but reads well 1. I have some problem with GDAL since, I think, I have installed OSGEO64 while I am using Python 2.7 32 bit. So I want to avoid GDAL.

Now what is the way to read 2 in arcpy.
I am using rasterio in virtual environment now and I want to get rid of the mess of using virtual environment.

My machine details:


texts: Network Analysis loop: fails to loop through desired points

texts:
This code is meant to automatically go through two sets of points, and determine distances from each point in one set to each point in the other set.  However,the code does not appear to actually cycle through the stations.  (When I test this, I should get 9 different distances, but instead each loop of three repeats the same set of numbers.  Also, the pStop point doesn't change when I test it.)

The code does work when i just link 1 station to 1 address.

stations =QgsMapLayerRegistry.instance().mapLayersByName('CTA_Stations_crs')[0]
addresses = QgsMapLayerRegistry.instance().mapLayersByName('Chicago Test Layer.shp')[0]

for train in range(0, 3):

    stationID = QgsFeatureRequest().setFilterFid(train)
    stationStation = next(stations.getFeatures(stationID))

    addressData=addresses.dataProvider()
    stationName = stationStation.attribute('LONGNAME')
    addressData.addAttributes([QgsField(stationName,QVariant.Double)])
    addresses.updateFields()

    for house in range(0,3):

        addressID = QgsFeatureRequest().setFilterFid(house)
        addressAddress = next(addresses.getFeatures(addressID))

        graphLayer = mmqgis_find_layer('Link.shp')
        director = QgsLineVectorLayerDirector(graphLayer, -1, '', '', '', 3)
        properter = QgsDistanceArcProperter()
        director.addProperter(properter)
        builder = QgsGraphBuilder(graphLayer.crs())

        pStart = QgsPoint(addressAddress.attribute('Longitude'), addressAddress.attribute('Latitude'))
        pStop = QgsPoint(stationStation.attribute('POINT_X'), stationStation.attribute('POINT_Y'))
        tiedPoints = director.makeGraph(builder, [pStart, pStop])
        graph = builder.graph()

        tStart = tiedPoints[0]
        tStop = tiedPoints[1]

        idStart = graph.findVertex(tStart)
        idStop = graph.findVertex(tStop)

        (tree, cost) = QgsGraphAnalyzer.dijkstra(graph, idStart, 0)
        Distance = cost[graph.findVertex(pStop)]/1840.18514097
        print(Distance)
        addresses.startEditing()
        attrIndex = addressData.fieldNameIndex(stationName)
        addresses.changeAttributeValue(addressAddress.id(), attrIndex, Distance)
        addresses.commitChanges()
        iface.vectorLayerTools().stopEditing(addresses)


texts: How to install ArcGIS 10.2 for Desktop on Windows 7 (64bit)?

texts:
I am going to install ArcGIS for Desktop 10.2 but I cannot find a tutorial about how to do that!

Does anyone have any advice on how to install this software?


texts: Leaflet : load & remove layer at specified interval time

texts:
I want to load and remove a layer in a specified interval of time. I make it blink with a random color function associate to my layer style. I can load this layer multiple times, unload it in a specified interval of time, but both actions at the same time failed (layer load and and load and load until Firefox crash). I also try clearLayer and setInterval(function() {getRandomColor();}, 200); without success.

How can I make my color function load in a specified interval of time ? (and make my style dynamic)

// Random colors
function getRandomColor() {
    var letters = '0123456789ABCDEF'.split('');
    var color = '#';
    for (var i = 0; i < 6; i++ ) {
        color += letters[Math.floor(Math.random() * 16)];
    }
return color;
}

// Layer style
function doStyleFields(feature) {
        return {
            color: '#000000',
            fillColor: getRandomColor(),
            weight: 0,
            dashArray: '',
            opacity: 1.0,
            fillOpacity: 0.5
        };
}

// Layer
function loadJson() {
    var exp_FieldsJSON = new L.geoJson(exp_Fields,{
        onEachFeature: pop_Fields,
        style: doStyleFields
    }).addTo(layerJson);
}

// LayerGroup
var layerJson = new L.LayerGroup();
var overlays = {"Fields": layerJson};
L.control.layers(overlays).addTo(map);

// Interval
$(document).ready(function() {
    setInterval(function() {loadJson();}, 200);
    setInterval(function() {map.removeLayer(layerJson);}, 199);
});


texts:
I am building a journey planner app using opentripplanner. i have followed their two minutes tutorial in which we just run webapp which uses pre build graph.obj file, it works fine.
But in next five minutes tutorial they shows how can i build .obj from a zip gtfs, while i try to do so i got these errors, i try had to figure them out !


texts: Exporting multiple shapefile layers from QGIS3 to one KMZ file

texts:
I have a .KMZ file which contains several separate layers which each contain important information, some are polygons, lines and points. When I import this into QGIS3 it displays separate KMZ layers organised into a group, even though it is only one single KMZ file.

I wish to edit some of these layers, and make new ones.

The issue is, when I come to export the KMZ after editing it seems like I can only export one layer at a time to KML. So I can open up a multi layered KMZ file in QGIS3, but I cannot export a multi layered (single file) KMZ file from QGIS3, I can only export single KMLs.

I wish to have one file which I can send to people to open up in Google Earth.

Is there a way to export multiple layers to KMZ in QGIS3? Or is there a work around where I can export each KML layer and then somehow combine them outside of QGIS3 to get me a KMZ layer?


texts: How to detect which Quickbird data product?

texts:
I have some Quickbird imagery, for further processing I would like to know what kind of data product are these images. The imagery only has an additional text file containing some order numbers.

Quickbird delivers Basic imagery were Orthorectification must be done. Ore they deliver Standard Imagery which is unsuitable for orthorectification.

It seems not to be a Ortho Ready Standard Imagery Product because with the data there was no RPC file.

Here is an example of [Readme.txt of Quickbird imagery][1]
Here is product description from Digital Globe.

[1]: ftp://96.53.39.210/pools/A/A0/WorldSatFTP/052761023120_01_README.TXT


texts: Workflow for sharing inventory of stormwater network in Open Source GIS?

texts:
Our agency has collaborated with six townships and one city for preparing an inventory of stormwater network.

We are using ArcGIS 10 (2 ArcView licenses) and ArcPad 7.1 (stuck with it since our Mobilemapper CX is running on Windows CE).

We are using ESRI File based Geodatabase (chose it because it allows multiple readers and one editor per feature datasets and it has one feature dataset for every political subdivision). We have finished 25% of the project.

My primary concerns are:

a) We cannot edit same feature dataset at a time. Mostly, my intern and I collect and edit data for same area/region at a time.

b) Collecting data in the field is not efficient enough since we cannot see each others edits. (using ArcPad Data Manager)

c) Sharing data with client is not a smooth task (converting the data into KML format so that it can be viewed on ArcGIS Explorer).

I understand that the solution to my problem is ArcSDE and ArcGIS server but we don't have funding. We want to finish the project ASAP.

Please recommend a workflow and open source software where we can work and share the data with clients more efficiently.


texts: How to get lines and nodes around the confluence point in a network system (line shapefile)?

texts:
I got a network line shapefile and the figure below was part of it. I want to do task around confluence point, so I took photo around confluence point as figure below.

I want to get line segments (unique geometry of line) and nodes (coordinates of node) around every confluence point in the shapefile, as the schematic figure below (closer look of figure above).

Now I use QGIS and manually select those segments and store node coordinates (via extract node) in CSV files. Because the network is too complicated and has many confluences, I want to do this in a much more efficient way.

The Idea
I think that maybe I could scan the whole network to find every confluence by finding node with more than 2 intersecting lines, and get set of line segments (in figure above, lines = line1, line2, line3) and nodes on the other side of the segment. (in the figure above, nodes = node1, node2, node3) (To be more specific, segment here means the line with detail in every row of attribute table of this line shapefile, not the little segment showed as I start editting.) Each confluence has its specific set, and in the case above, the set is:

geom of line1, coordinates of node1
geom of line2, coordinates of node2
geom of line3, coordinates of node3

How can I get the line segment and its node around the confluence point (auto-determine-and-find) in a network-like line shapefile?

And as I described in the paragraph above, because line shapefile I got is complicated, there could be many confluence points, thus I want found lines and nodes could be grouped in the same set based on the confluence points which they intersected.

I don't know if my idea works and what function or tool could help me to achieve my goal. I want to achieve this via open-sources, and ultimately write a integrated python function to get those lines and nodes.


texts: Enable Popups for each layer in webmap using ArcGIS API for Python

texts:
I am looking to use the ArcGIS API for Python to enable popups for each layer in a webmap.
I found Configure Popup Attributes Programmatically with ArcGIS API for Python on GeoNet which addresses how to disable popups for each layer, but I am looking to do the opposite.
I have also asked Enable Popups for each layer in a webmap Enable Popups for each layer in a webmap using ArcGIS API for Python on GeoNet.


texts: infer the altitude from a DTM & finally provide it to point / evincere la altitudine da un DTM & finalmente fornirla a dei punti

texts:
I would like to complete a point vector's table with Z data.
I am not able to take Z coordinate from a dtm
(downloaded from  http://www.regione.piemonte.it/geoscaweb/index.php?pgs2=0&pgs3=DTM10&pgs4=DTM10_f&pgs5=.zip&pgs6=%20DTM%20(passo%2010%20m)%20-%20Modello%20altezze%20filtrato%20(fogli%2050.000)&pgs8=50),
   I'm using GRASS from QGIS. All files are in its database, correctly imported.

I'm trying with v.what.rast but , even if it tells he succeeds, no new data has been recorded in the point vector's table and no new column appears.

I'm questiong about
1)the possibility of doing that with this particular DTMfile;
2)the fact that there aren't specific requirements for the point vector's table.


texts: Calculating walking distances on network in ArcGIS Network Analyst?

texts:
I would like to calculate distances on network between demand nodes (points) and transportation network (line) in ArcGIS.

I basically calculated the central features of my building data (1,5M feaures) in considered district zones (85 zones). Now i have 85 central node. I also the highway network.

I know i can calculate the distance if i use the euclidean distance tool. However, for representing the real walking distance, i need to find the real walking distances on the network.

What i want to do is calculate these demand node points' rectilinear distances on network instead of euclidean distances. I need only one distance value for a demand center to a highway line.

What should i do?


texts: Error "std::exception" in GRASS GIS via QGIS

texts: Does anyone know what this error means? I am using the GRASS GIS plugin in QGIS and I have imported a road layer (vector) using v.in.ogr.


texts: Exchange Boundaries and Wire Centers in France and Germany

texts:
I am attempting to find a shapefile of wire center service areas / telephone exchanges in France and Germany.  Ideally, a point file showing the location of the actual equipment, coupled with a polygon that shows each center's service area.

To see an example of the what I'm talking about, this is the U.S. equivalent:
http://www.navmart.com/wire_center_boundaries.php

Obviously, free is best, but I open to sources who might have this available for purchase.

Any and all help is appreciated!


texts: Adding CSV file from Google Drive to ArcGIS online Map Viewer

texts:
I am trying to add a layer from web in ArcGIS Online's Map Viewer. I had uploaded a csv with proper location details (Long , Lat) in my Google Drive and made it public shareable link.
When I try to import it to ArcGIS Online it says

The layer was not created because no locations could be found. Make
sure your data has valid location information.

Is there a solution for this?


texts: I would like to know is there anyway to overline a text label like this in ArcGIS10.3.


texts: Need to change projection and data to match 2nd shapefile

texts:
I have two files I am using, one is census tract data with prj:

GEOGCS["GCS_North_American_1983",
    DATUM["D_North_American_1983",
         SPHEROID["GRS_1980",6378137.0,298.257222101]],
    PRIMEM["Greenwich",0.0],
    UNIT["Degree",0.0174532925199433]]

The other is a grid file from a model with prj:

PROJCS["airpact",
    GEOGCS["GCS_Sphere_EMEP",
       DATUM["D_Sphere_EMEP",
         SPHEROID["Sphere_EMEP",6370000.0,0.0]],
       PRIMEM["Greenwich",0.0],
       UNIT["Degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic"],
    PARAMETER["False_Easting",0.0],
    PARAMETER["False_Northing",0.0],
    PARAMETER["Central_Meridian",-121.0],
    PARAMETER["Standard_Parallel_1",30.0],
    PARAMETER["Standard_Parallel_2",60.0],
    PARAMETER["Scale_Factor",1.0],
    PARAMETER["Latitude_Of_Origin",49.0],
    UNIT["Meter",1.0]]

Data file
In layer spatial reference system units

xMin,yMin -124.849,41.988 : xMax,yMax -111.043,49.0025

Layer Spatial Reference System

+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs

Grid file
In layer spatial reference system units

xMin,yMin -342000.00,-942000.00 : xMax,yMax 798000.00,90000

Layer Spatial Reference System

+proj=lcc +lat_1=30 +lat_2=60 +lat_0=49 +lon_0=-121 +x_0=0 +y_0=0 +a=6370000 +b=6370000 +units=m +no_defs

My problem is obvious, these are nowhere near the same and when I load into qgis the grid is far larger than the census tract data. I need to define them both on the same projection as the grid without ruining the census data but I am not sure how to do this. I have read several similar questions but am still lost.

I am aware that the census data is referencing lat/lon in degrees and the grid is in projected coordinates in meters, but need to somehow convert the census data to match the grid.


texts:
Issue 1: I have been working on the re-evaluation of a Late Woodland Palisaded Village in southwest VA. The site map as been digitized in CAD. The original grid has been followed. The excavators did not start their grid in the NE or +,+ quadrant, but instead the grid is in the NW or -, + quadrant. All of the features have been traced using a CAD polygon tool and all have been numbered with an unique number. So far so good.

BUT, as it turns out, the excavators had three levels with vital information apart from the feature shapes. So all 1011 ten foot squares have been digitized using a duplicate array tool. And each composite level set can be copied and pasted exactly into place using standard CAD tools. At which point, manual labeling of each and every square of the two other level would be required.

The problem is how to assign labels according to the oddball grid system such that when the 3033 squares are shapefiled for import into QGIS, one can do useful things with those squares.

Being a relative newbie to GIS, it struck me that the great divide between traditional archaeology, of which perhaps a majority of practitioners use CAD; and GIS, where each square is geo-located, is in how to label each square such that when turned into shapefiles, the labels are in the attribute tables and can be joined to other data. The Z values are not important as I am trying to show distributional data and then can merge them after searches. My CAD solution is for each of the first 1011 squares, I assign an alpha-numeric code to each that corresponds to the original map. The entertaining factor is that the traditional X,Y system is reversed in that the Y or column is listed first, then the X or row. The X values are numeric and the Y values are Alpha. Either way, my problem is to either plow ahead spending days for all of these squares, or whether there is a far simpler and somewhat automated macro/scripting system to accomplish that?

So, can GIS, any of the flavors, although I use QGIS, accomplish the automated meta-data assigning of an arse-backwards system of recording to a set of squares? To further complicate matters, not every square was excavated. But, if there were a system that allowed for a rectangular array of squares to be coded, and then imported into CAD, I could easily enough remove the not needed squares by manual selection and deletion.

Rounding back onto the original problem, it appears that getting flat maps into geo-coded GIS applicable systems is the biggest hurdle that combining archaeology and GIS faces. If there were a relatively easy way of geo-coding features/squares/whatever such that they could then be imported into GIS, a lot more useful analysis might be performed.

But, as a relative newbie, am I totally off-base and is there a relatively easy solution?


texts:
We are using QGIS 3.4.14 -Madeira - Vector Tiles Reader plugin (3.2.2). What version of MapBox GL style JS is used in the Vector Tiles Reader plugin?
Regards,
Rajesh


texts: How to Prevent Final (Unnecessary) Selection Box in Script Tool Input Parameters via Validation Script?

texts:
Using Arc 10.2.2.

My script tool's interface has its first parameter populated by the following snippet of validation script:

  def initializeParameters(self):
    # Source Folder
    dir_f = r'Q:\GW\\EC1210WQAEH_QESEA\\CSSP_ATL\\Data\\CSSP_Web_Tools_Imports'
    csvs = [os.path.join(dir_f, f) for f in os.listdir(dir_f) if f.endswith('.csv')]
    self.params[0].filter.list = csvs

It populates the first parameter window with all the .csv files in the specified directory but for some reason it adds a final selection box and makes it checked.

Screen grab:

Every time I run the tool I must scroll to the bottom of the list and uncheck this box to make the red "X" (error message) go away and allow me to run the tool.

Why is this happening?


texts: Can multiprocessing with arcpy be run in a script tool?

texts:
I have been tinkering around with the multiprocessing module to resolve some processing of data in ArcMap 10.3.

I have created a script which works fine if I run it in IDLE. I see all my cores max out in Task Manager and the code completes without error.

Now if I wire this script up as a Script Tool in ArcToolbox it throws a weird error

Could not find file: from multiprocessing.forking import main; main().mxd

Now reading various threads I regularly see that the script has to be "run out of process". So I unticked the check box "run python in process" in the script properties dialog and was expecting the code to run but it does not, I get a 000714 error and I've done nothing to the code.

So my question is simply can one create a script that uses the multiprocessing module and run it from ArcToolbox? From my limited playing around, it seems I cannot and its a technique that can only be run from IDLE?

I have finally got it working, using every ones help. I decided to fully document my simple example and bring together any of the trip ups I had to overcome in a document on the ESRI GeoNet website. Like Luke's answer I hope this provides a starting point for anyone trying to create a python script tool that utilises multiprocessing. The document is titled Create a script tool that uses multiprocessing.


texts:
I have a full OSGeo4W (64bit) installation on my Windows 10 laptop.
I'm using it for Quantum GIS, Python et cetera. There is a osgeo4w64/apps/gdal-dev folder in the installation.

I want to achieve this; http://www.cruisersforum.com/forums/f134/opencpn-13653-50.html#post322600

Which installs some extra packages ('frmts' from here; https://github.com/OSGeo/gdal/tree/trunk/gdal/frmts/bsb) and makes use of it. Now my question(s);

what are these 'frmts', some sort of plugins for GDAL?
Are they already installed on my OSGeo distribution?
If yes, where, so I can uncomment the line.
How to download and install this package?


texts: Calculate elevation from rgb from png created from qgis

texts: I am using qgis to create a raster elevation image (png) from SRTM data. I need to create a png via qgis which I will then query to determine the rgb of particular pixels via php. How do I define a color scheme within qgis in a way that I can then use the same gradient to calculate the elevation of the pixel from the png?


texts: I have a DEM that I'd like to render in 3D using QGIS 3.12. Various tutorials tell me to click on view>3D map view to do this. However, when I do this no window or anything pops up. If I go to View>Panels, it says that I have a panel for a 3D map view open, but it does not appear on my screen. Have tried updating my graphics drivers and reinstalling QGIS.


texts: Is it possible to disable infowindow but retain featureClick event in cartodb.js?

texts:
I would like to pass the data from a marker in a featureClick() event to perform other actions. However I do NOT want an infowindow or tooltip to display. I can quite easily remove the infowindow/tooltip in the options of createLayer(), however this disables the click functionality for a marker.

EDIT 1:
It seems that if you explicitly setInteraction(true) for the sublayer the featureClick() event will fire, however the cursor for the marker is not changing to a pointer.


texts: How do I apply Dijkstra's algorithm in Python with PostGIS data?

texts:
I have currently parsed OSM data, from a PBF file, into a pgRouting database with the use of osm2po, an excellent parsing tool. I don't however want to use the algorithms provided by pgRouting and instead want to have a custom algorithm. But to start off and gain an understanding of how to use the data, I'll attempt to simply apply Dijkstra's algorithm.

My logic so far is to:

Determine whether two points are within x distance.
Get the two nodes closest to each point.
Query the database with a bounding box.

The query to the database will return all nodes, including their: x and y position, target nodes, edge length (km) and the_geom, within the specified bounding box.

The difficulty I have now is, how can I use the data and apply Dijkstra's algorithm such as the algorithm here: http://thomas.pelletier.im/2010/02/dijkstras-algorithm-python-implementation/, as the query won't produce a graph like:

...      graph = {
...     'A': {'B': 10, 'D': 4, 'F': 10},
...     'B': {'E': 5, 'J': 10, 'I': 17},
...     'C': {'A': 4, 'D': 10, 'E': 16},
...     'D': {'F': 12, 'G': 21},
...     'E': {'G': 4},
...     'F': {'H': 3},
...     'G': {'J': 3},
...     'H': {'G': 3, 'J': 5},
...     'I': {},
...     'J': {'I': 8},
...     }

If Dijkstra's algorithm is not suitable for this, what routing algorithm is?

In addition to that, I'm slightly confused how routing algorithms work on parsed data, which has been parsed by any osm parser. From all the parsers I've come across, a node only ever has one target node. Wouldn't this mean that I could only ever move in one direction, that being to the target node, then the next target node and so on?

Thanks to anyone willing to improve my logic.


texts: Is it better to define object properties or layer properties?

texts:
I have a question, but it is difficult to phrase. I am picking up GIS methods as I go, and I have run into this problem frequently. Please excuse me if I am in the wrong place.

When organizing information, is it better to group Features together and display/arrange accordingly? Or is it better to attach specific information to each Feature individually?

For a real-world example, I have lots of Features (>10000). Each Feature is given an icon, where the icon represents aspects of the feature. The map consists of layers that organize and display the Features. Take, for example, the question of where to store the icon information? Should I let the layers dictate the icons that are displayed? Or should I let the Features render themselves in whatever layers they are in? This is coming from a programmer, so processing speed/memory are the main issues.

Is there a general design principle for this sort of behavior? If my question indicates a deeper misunderstanding, please fill me in. I've been running in circles.


texts: Active Directory Authentication & MS SQL Server in ArcReader

texts:
I'm having difficulty with a single user (That I know of) being able to access geodatabase files in our ArcGIS Server using ArcReader.

Our setup is as follows: Enterprise Geodatabase within MS SQL Server 2012, Active Directory integration within SQL Server (GIS Edit group & GIS View group), mixed environment of ArcGIS users and viewers.

Every other GIS user with either ArcReader or ArcView can access the geodatabase data, but this new user I have added to the ArcView group cannot see the layers. The PMF map shows all of them as unavailable. Her ArcReader has the server connection in the proper folder to allow for OS authentication, and I even tried putting her in the edit group but she still cannot access the layers.

We're running 10.2 server and she has 10.2.1 ArcReader.


texts: ArcGIS for Android API Example/Tutorials

texts: I know that the ArcGIS API for Android is still in beta (I've joined the beta program), but I was curious if anyone had produced or discovered any examples/tutorials/cookbooks for the API. The Help Resource is still severely limited. I have watched the ArcGIS API for Android Videos, but they do not provide any of the source code (contrary to what is stated in the videos themselves). I have also browsed through the API Reference but I've found the documentation to be lacking.


texts: How to Specify Haversine when using Buffer Method in Shapely and how to get Haversine distance between two Shapely Point objects?

texts:
I have a series of latitude/longitude pairs and I form Point() objects for each pair.  I would like to form buffer regions around each point using a Haversine distance argument instead of Euclidean distance.  That way, points on the boundary of the buffer are an equal Haversine distance away from the center of the buffer.  Is there a way I can do this?

Also, is there a built-in method for getting the Haversine distance between two Shapely Point() objects?  I know I can use the Haversine formula for the lat/long pairs, but I was wondering if there was anything built-in.


texts:
Following the advice on Installing QGIS with GRASS plugin on Ubuntu?...

I've installed qgis, qgis-plugin-grass, python-qgis, grass-core=6.4.3-3, grass=6.4.3-3, grass-doc=6.4.3-3, grass-gui=6.4.3-3

I open qgis, but there is nothing about GRASS in the plugins menu.

scott@beast:~$ apt-cache policy qgis
qgis:
  Installed: 2.8.3-0~trusty1
  Candidate: 2.8.3-0~trusty1
  Version table:
 *** 2.8.3-0~trusty1 0
        500 http://ppa.launchpad.net/ubuntugis/ppa/ubuntu/ trusty/main amd64 Packages
        100 /var/lib/dpkg/status
     2.0.1-2build2 0
        500 http://us.archive.ubuntu.com/ubuntu/ trusty/universe amd64 Packages
scott@beast:~$ apt-cache policy grass
grass:
  Installed: 6.4.3-3
  Candidate: 7.0.1-1~exp1~trusty
  Version table:
     7.0.1-1~exp1~trusty 0
        500 http://ppa.launchpad.net/ubuntugis/ppa/ubuntu/ trusty/main amd64 Packages
 *** 6.4.3-3 0
        500 http://us.archive.ubuntu.com/ubuntu/ trusty/universe amd64 Packages
        100 /var/lib/dpkg/status


texts: SLD : setting Multiple Values in one Literal Tag

texts:
Is it possible to use multiple values in one Literal tag?

This is what I have:

    <ogc:PropertyIsEqualTo>
      <ogc:PropertyName>name</ogc:PropertyName>
      <ogc:Literal>0</ogc:Literal>
    </ogc:PropertyIsEqualTo>

I want to symbolize multiple attributes in one block rather than repeat the above line dozens of times.  And unfortunately my table values are not very friendly so I can't use a greater than call.  I was hoping for something similar to below (Literal Tag)

    <ogc:PropertyIsEqualTo>
      <ogc:PropertyName>name</ogc:PropertyName>
      <ogc:Literal>0;10;60;40</ogc:Literal>
    </ogc:PropertyIsEqualTo>

And the next Rule would be something like this:

    <ogc:PropertyIsEqualTo>
      <ogc:PropertyName>name</ogc:PropertyName>
      <ogc:Literal>1;11;62;42</ogc:Literal>
    </ogc:PropertyIsEqualTo>


texts: Rule-based labelling in QGIS - problem moving labels

texts: I have just started to use the Rule-based labeling feature in QGIS but am not able to move the labels - x and y coordinates are defined


texts: Populating values in field based off of another sorted field?

texts:
I'm working with a dataset of center pivot irrigation sites for my work and I have a unique ID for each land owner (Land_Owner_ID) as well as a number of center pivots associated with each owner (CP_ID). Ignore the certificate number, as it is irrelevant to my issue.

For each unique Land_Owner_ID I need to reset my CP_ID to 1 for the 1st instance of the unique Land_Owner_ID, then add 1 to CP_ID and set that as the next CP_ID for the same instance of Land_Owner_ID.

As a result, in the case of Land_Owner_ID 3 I should have CP_ID's 1, 2, and 3, and then again for Land_Owner_ID 4 I should have CP_ID's 1, 2, and 3.

I am somewhat of a python beginner, but I've been trying to do as much of my work via python scripts/tools as possible to flesh out my knowledge and experience.

Here is what I have come up with so far:

#This script sorts the field LandOwner_ID low to high and then populates the
# CP_ID and CP_LandOwner_ID_Combined fields
import arcpy
import os
from arcpy import env
env.overwriteOutput = True
arcpy.env.workspace = r"D:/Irrigation_Data/WaterUseData.gdb"

CP_shp = "All_CP_With_Masterlist_Version1_Copy"
out_path = arcpy.env.workspace

out_shp = os.path.join(out_path,"
 {}".format("All_CP_With_Masterlist_Version1_SortedCopy"))
sort_field = [["Land_Owner_ID","ASCENDING"]]

#Sort Land_Owner_ID from low to high
arcpy.Sort_management(CP_shp,out_shp,sort_field)

landOwner_ID_Field = "Land_Owner_ID"
CP_ID_Field = "CP_ID"
updateFields = (landOwner_ID_Field,CP_ID_Field)

with arcpy.da.UpdateCursor(out_shp,updateFields) as cursor:
  for row in cursor:
     #Reset the CP_ID for each new landOwner_ID to 1
     CP_ID_Count = 1
     #Store owner_ID in a variable
     landOwner_ID_Storage = 1
     landOwner_ID = row[0]

If current landOwnerID equals the ID I am testing for, while that condition is true, set the current CP ID and count up 1, continue to do so until a new landOwner_ID is registered. I'm sure I have some syntax or arrangement of my code out of whack here.

     if landOwner_ID == landOwner_ID_Storage:
         while landOwner_ID == landOwner_ID_Storage:
             row[1] = CP_ID_Count
             CP_ID_Count += 1
             cursor.updateRow(row)

     else:
         landOwner_ID_Storage += 1


texts: Measure for spatial autocorrelation in fire occurrence data

texts:
I am trying to determine the variables that predict forest fire occurrence in a region in Chile.

Firstly, I determined burned areas using LANDSAT data and the dNBR-index. I ultimately made a raster with 30m forest pixels with a value of 1 (burned) or 0 (not burned).
I wanted to do a random sampling with minimum distance to reduce the autocorrelation between burned pixels.
My problem is that I do not know to what distance my burned LANDSAT pixels are autocorrelated.
I tried to make a semivariogram in R but this gave errors, probably due to the binary variable or because there are pixels with no-data values. I made a point shapefile of my raster for which every point has the same binary value of the pixel it is covering. However, if I make a semivariogram, I assume that the autocorrelation of pixels that were not burned isconsidered as well and this is not what I want.

Anyone an idea to investigate the autocorrelation of points or pixels with a specified value (1 in my case, which represents fire).


texts: Popup - a various set of attributes depends of selected layer

texts:
I would like display in popup a various set of attributes depends of selected layer.
I have:

    var selectControl =  new OpenLayers.Control.SelectFeature(
    [layer1, layer2],
    {
    clickout: true, toggle: false,
    multiple: false, hover: false,
    toggleKey: "ctrlKey", // ctrl key removes from selection
    multipleKey: "shiftKey" // shift key adds to selection
    }
);

and:

    function onFeatureSelect(event) {
    var feature = event.feature;
    var selectedFeature = feature;
    var text = ""
        if(feature.cluster) {
            for(var i=0; i<feature.cluster.length; i++) {
//and here should be something about which layer is selected eg.
//if ( selectControl == layer1 )
//then {first set of attributes} else {second set of attributes}
//... How to do htis?
//This idea doesn't work, I received undefined value of attributes

text += [i+1]+". "+feature.cluster[i].attributes.gis_id +"<br>";
}
else {
    text = feature.attributes.gis_id;
    }
    var popup = new OpenLayers.Popup.FramedCloud("chicken",
        feature.geometry.getBounds().getCenterLonLat(),
        new OpenLayers.Size(100,100),
        text,
        null, true,
        onPopupClose,
        function() { controls['selectControl'].unselectAll(); }//
    );
    feature.popup = popup;
    feature.popup.closeOnMove = true;//
    map.addPopup(popup);
};


texts: Anyone knows why is this happening? It does not let me to save my script.


texts: Writing 3 channels to 8-bit TIF in Python using gdal

texts:
I have some Python code to write 3 channels to file. Here, channels is an array of 3 numpy 2D-arrays (RGB), outfile is the filename, and rows and cols are the image dimensions.

def save_tiff(channels, outfile, rows, cols):
    outdriver = gdal.GetDriverByName("GTiff")
    outdata = outdriver.Create(str(outfile)+".tif", rows, cols, len(channels), gdal.GDT_Byte)

    # write the arrays to the file
    i = 1
    for c in channels:
        outdata.GetRasterBand(i).WriteArray(c)
        outdata.GetRasterBand(i).FlushCache()
        i += 1

However, this results in 24-bit output images (probably because of 3xGDT_Byte channel). How do I get a single 8-bit image from 3 channels (R,G,B) in Python using gdal?


texts: Making square footprint polygons in QGIS

texts:
We are a group of students working with footprint taken from LiDAR point cloud. We have a challenge in making the polygons more square. We want the angles to be orthogonal, and the sides to be longer than eg 3 meters. We have tried QGIS "generalized" and "simplified", but are not completely satisfied.
Does anyone know of a good tool for this task?


texts: Reclassify a raster based on proximity to other cells in the same raster

texts:
Using Arc 10.2.2 with the Spatial Analyst extension, I have a single band, signed integer, 8-bit raster.  I need to reclassify cells within that raster based on their Euclidean distance to other cells within that same raster.

For example, if a cell with value 2 is within 100 meters of any cell containing value 5, then reclassify the value 2 cell to value 3.  And so on.

I considered FocalStatistics, but none of that tool's statistic types (mean, majority, etc) met my needs.

A search on the topic turns up examples where the proximity is to features within a separate vector layer, but nothing on proximity within the same raster.


texts: Creating fishnet grid Shapefile in QGIS?

texts:
How do I create a fishnet shapefile using QGIS?

The equivalent task of "Create Fishnet" in ArcGIS.


texts: Shapefiles of Central America and Mexico coastlines?

texts:
Where can I find a shapefile of Central America + Mexico (WSG84) (with no administrative boundaries) at a low resolution?

I'm looking for a very simple shapefile without the small islands and smoothed out coastlines.


texts: Clip polygon features with country boundary ArcGIS Desktop

texts: I am trying to trim the edges of some polygons based on a country boundary. I have tried both clip and intersect, but neither seem to be doing what I want. Below is the image I am working with. When I use clip and intersect, I only get the purple polygon returned. Should it not trim the edges of the other polygons to inside the country boundary, or am I using the wrong tool?


texts: Materialized View Does Not Appear In Geometry_Columns View

texts:
I have defined a materialized view in a PostGIS 2.1-enabled database running PostgreSQL 9.3.1:

CREATE MATERIALIZED VIEW canvec_data.contours_m AS
 SELECT
    contours_original.ogc_fid,
    contours_original.elevation::integer AS elevation_orig,
    CASE
        WHEN "substring"(contours_original.code::text, 1, 3) = '257'::text THEN round(contours_original.elevation::double precision * 0.3048::double precision)
        ELSE contours_original.elevation::double precision
    END AS elevation_m,
    contours_original.wkb_geometry::geometry(Linestring,3578)
 FROM canvec_data.contours_original

I had, perhaps erroneously, expected that this resultant entity (a table?) would be listed in the public.geometry_columns view, and thus be properly interpreted by GIS viewer software such as QGIS.

The materialized view does not appear in the public.geometry_columns view.

What am I doing wrong?

Thanks!


texts: Why Haversine formula and Web Mercator give different length between two latlng points?

texts:
I am calculating the distance between two points given by latlng. However, Google Maps (i.e., web mercator projection) gives me a longer distance than Haversine formula.

Which is the most accurate measure? I can read in some answer that Haversine formula gives it.

So... Is Google Map lying when telling me distance in meters. e.g., in this webpage? http://www.maptiler.org/google-maps-coordinates-tile-bounds-projection/

Edit:

Example: In a small area around (37,0), if I calculate area (width,height), haversine results in (7750m*7753m) while Google results in (9784m*9784m).


texts: ArcSDE reconcile versions and keep edit version in its current state

texts:
Reconciling versions makes both target and source(s) point to the same state.  I have an arcpy script that essentially will create a version from a parent of the db - in a state that is empty - adds data to the feature classes, and then will append them to its parent's parent.

So its immediate parent is still pointing at a state that is empty.  But now the child and the one it is merging with(state above its parent) are both point to a version that has a lot of data.  Ideally, I would like a way to keep the child version just have the data that was loaded to it.  Essentially have the child contain only data that came in during that process.

This seems like it should be a thing and would be a quick way to copy a lot of data to a version without the use of any cursors.  I found this post
When Versioning with ArcSDE can posted edits be cancelled or rejected?
and it talks about editing the sde versions table.  Hoping there is another option.  However is there a clean way of changing the state using ArcPy? I am thinking, if I have to, I will just connect with pyodbc and then change the table.

Also, is it 100% that when a database version is created and then immediately reconciled without any other users or transactions happening, that its state will only one ahead of its creation?


texts: Why measure height relative to geoid instead of ellipsoid?

texts:
Why do we measure elevation relative to sea level? I understand that water settles to a level surface, but why not just measure relative to an ellipsoid and dodge all the sea level rising, geoid, and other problems. I'm looking for a real world example where height relative to a geoid is more useful than height relative to an ellipsoid. What could go wrong if I used height relative to WGS84 ellipsoid for example?
I'm new to all this and thought about it for a day but couldn't figure it out.


texts:
I need the bounds/full extent of a WMS layer in Leaflet. I don't think L.TileLayer.WMS has a getBounds() method. I had this issue with L.ImageOverlay, and added the functionality via:

L.ImageOverlay.include({
    getBounds: function () {
        return this._bounds;
    }
});

Is there a solution for L.TileLayer.WMS? Or will I have to start parsing the WMS getCapabilities request for its LatLongBoundingBox or BoundingBox elements?

EDIT

I've just found a duplicate: Get bounding box of WMS. This question is still relevant, since Leaflet 1.0 has since been released.

EDIT 2

Get bounding box of WMS does not answer my question - its answer is to parse the getCapabilities response manually. Is there any other way?


texts:
I use following script to determine text size of different labels:

CASE
WHEN "KKOD" IN (1) THEN 8
WHEN "KKOD" IN (2) THEN 9
-- etc...(here come more lines of similar code)

ELSE 0

END

Everything works fine. But now I want to scale all attributes to for example 150%. Is there any simple way to do it?
Is it possible to put attributes in some variable and then multiply it by 1.5 (150%)? Something like this:

"KKOD" * 1.5

I don't want to change all numbers manually because it will be difficult to change them to something else in the future.


texts: How to fix "The entry point for the procedure sqlite3_open_v2 cannot be found in the dynamic link library sqlite3.dll"

texts:
Translation:
Hello, I have a problem installing version 1.0.8. The following error message appears:

The entry point for the procedure sqlite3_open_v2 cannot be found in the dynamic link library sqlite3.dll

[The next sentence is ungrammatical in the original and therefore its translation is in doubt.]
Previously I had version 1.7 Wroclaw, which I uninstalled before installing the new one.
Thank you for your help.
Original
Bonjour,
J'ai un problème pour installer la version 1.0.8, le message d'erreur suivant apparait:

Le point d'entrée de procédure sqlite3_open_v2 est introuvable dans la bibliothèque de liens dynamiques sqlite3.dll

J'avais avant d'installer la version 1.7 Wroclaw que j'ai désinstaller pour installer la nouvelle.
Merci pour votre aide.


texts: GRASS PyWPS Create a new vector map from user selected polygon

texts:
My goal is to select a specific polygon in which the user clicked and make a new vector map that contains only the selected polygon.
I see two ways:

Convert XY coordinate from literal input to a vector point map and use v.select. But I can not find a command (or a workaround) for the first part?
Use v.what (XY from literal input) to see which polygon is located at this point and use ID information in v.extract command as binded variable with SQL where query. But how to store/use for example, "ID" data (printed in shell script style output from v.what) as global variable?

Also, the idea is first to use standard GRASS modules with cmd.