texts:
I'm trying to create a shapefile I can use as a template for a breeding bird survey. The basic gist was that each species has a two letter code and a colour depending on the BTO listing. Each bird is either red, amber or green. There would then be extra markers to say whether the bird was singing or calling etc. So, I set up an attribute table with a column labelled "BTO Code" and "Behaviour"

I then went ahead and set it up so that when you typed a species code into that column, it would set the point marker as those letters coloured either green, amber or red. At first it was working, I kept testing it as I went along and the points were showing up fine on the map. The markers for singing/calling were working fine too. Below is what shapefile and my rules looked like:

Well, I just finished adding the last few species, clicked apply and ALL my points disappeared. Now when I try to add a species, they appear in the attribute table but not the map!

I'm at a loss. This is my second attempt at creating this shapefile/styling and the same thing happened last time, last time I assumed it was something I had done wrong, and couldn't for the life of me figure out what it was so I just started again. I really don't have the will to do this all over again! I'm at my wits end. Please can someone with a bit more knowledge than me shed some light on why it would suddenly stop working when I've done everything right in the code?

All I can think is that there may be a limit to how many rules you can set when using rule-based styling? There are a lot of birds! I've looked online though and can't find anything regarding a limit to how many rules you can set.


texts: Export as DGN with QgsVectorFileWriter.writeAsVectorFormat

texts:
This code works in QGIS 2.14 but not in 2.18. What do I have to change to be able to export properly to DGN. From what I understand, now it tries to save in microstation V8 format instead of v7.

I don't want to use "convert format" from the processing tool box, I'd really like to save directly to DGN without intermediary steps.

QgsVectorFileWriter.writeAsVectorFormat(
    couche_source,   #loaded layer
    enregistre_sous, #path to save as
    'utf8',
    QgsCoordinateReferenceSystem(2950, QgsCoordinateReferenceSystem.EpsgCrsId),
    'DGN',
    False,
    None,
    skipAttributeCreation = True,
    )

I really don't know what to do, and I didn't find anything helpful to me.

Even when I try to save as DGN directly form QGIS it doesn't work. I get this error message :

Export to vector file failed.
Error: creation of layer failed (OGR error:creation option '' is not formatted with the key=value format)

And I did not modify anything in the options.


texts: Differences between APP-6A and APP-6B in LuciadMap?

texts:
I'm using luciadmap 2013.1 and, from what I understood, one can switch between APP-6A and APP-6B symbology like this:

new TLcdEditableAPP6AObject(ELcdAPP6Standard.APP_6A)
new TLcdEditableAPP6AObject(ELcdAPP6Standard.APP_6B)

Now, I was wondering about the differences between both specifications and I wanted to display a symbol that would look different in APP-6A and APP-6B.

Does anybody know the code of such a symbol? And if it does not exist, then maybe the code of a symbol that appeared in APP-6B.


texts: Measure feature length of a shapefile that falls within boundaries of polygon layers of a different shapefile

texts:
I use Qgis 2.18.6

I have a vector layer with roads. These are of two types, new and existing roads, in the same shapefile.Those roads, on the map cross over different propertys (propertys of land).

I have a raster layer from a different shapefile with those property lines.

I want to know how long the roads are on each of the propertys. Preferably as a table.


texts:
I have a problem with the grass function r.mfilter outputting an empty file. I'm running the tool from qgis processing standalone.

My system:

OS: Ubuntu 16.04.2 LTS
Saga: 2.2.3
Grass: 7.0.3
QGIS: 2.18.4

Saga installed as per instructions in this question

The script:

from qgis.core import *
from qgis.gui import *
import os
import sys
from PyQt4.QtGui import *
from PyQt4.QtCore import Qt

app = QgsApplication(sys.argv, True)
QgsApplication.initQgis();
QgsApplication.setPrefixPath('/usr', True)

sys.path.append('/usr/share/qgis/python/plugins')
from processing.core.Processing import Processing
import processing.tools.general as processing

Processing.initialize()

# Processing directory
cwd = os.getcwd()
procdir = cwd + '/LST-trad/hojd/'
tmpdir = cwd + '/LST-trad/Test2/'
outdir = cwd + '/LST-trad/treecover2/'
filter = cwd + "/LST-trad/qgis5x5filter_1or0_div1.txt"
files = os.listdir(procdir)
# Filter files
filefilt =[]
for f in files:
    if f.endswith(".tif"):
        filefilt.append(f)

print "help:"
#processing.alghelp("saga:rastercalculator")
processing.alghelp("grass:r.mfilter")

# Coordinates
arglist = []
for file in filefilt:
    x1 = int(file[7] + file[11:13]+"000")
    x2 = x1 + 25000
    y1 = int(file[4:6] + file[9:11] +"000")
    y2 = y1 + 25000
    coords =  str(x1) + "," + str(x2) + "," + str( y1) + "," + str(y2)
    outfile = "COV_" + file[4:]
    arglist.append([procdir+file, coords, outdir+outfile, outfile])

#print arglist

formula1 =('ifelse(a > 50, 1, 0)')
formula2 = ('a * 4')

for arg in arglist:
    tmp1 = tmpdir + "t1_" + arg[3]
    tmp2 = tmpdir + "t2_" + arg[3]
    print "saga"
    processing.runalg("saga:rastercalculator", arg[0], None, formula1, True, 1, tmp1)
    print "grass"
    processing.runalg("grass:r.mfilter", tmp1, filter,1,False,arg[1], 0, tmp2)
    print "rascal"
    processing.runalg("saga:rastercalculator", tmp2, None, formula2, True, 1, arg[2])

QgsApplication.exitQgis();

The filter file:

MATRIX    5
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
DIVISOR   1
TYPE      P

output form alghelp(r.mfilter):

ALGORITHM: r.mfilter - Performs raster map matrix filter.
    input <ParameterRaster>
    filter <ParameterFile>
    repeat <ParameterNumber>
    -z <ParameterBoolean>
    GRASS_REGION_PARAMETER <ParameterExtent>
    GRASS_REGION_CELLSIZE_PARAMETER <ParameterNumber>
    output <OutputRaster>

One thing that confuses me is that the parameters listed here does not match the parameters listed in the grass documentation


texts: Enable/Disable Parameters based on Boolean

texts: I'm working in model builder and I have a boolean parameter that determines which of two processes run.  Right now, all of the parameters for both processes show up when I run the model.  Is there a way to only enable the parameters that are necessary for the one process that will run based on the boolean?


texts:
I am having difficulty understanding the following error message.

    var okombahe_erven = L.tileLayer.betterWms("http://localhost:8080/geoserver/Erongo/wms", {
    layers: 'okombahe_erven',
    format: 'image/png',
    transparent: true,
    attribution: '',
  });

  var map = L.map('map', {
    zoomControl:true,
    zoom: 14,
    maxZoom:25,
    minZoom:5,
    center: new L.LatLng(-21.3665370941162, 15.3731527328491),

    layers: [OSM, okombahe_erven]
  });

  var baseMaps = {
    "OpenStreetMap": OSM,
    "ESRI  World Imnagery": Esri_WorldImagery
  };
  var overlayMaps = {

    "Zonning": okombahe_erven,
  };
  L.control.layers(baseMaps, overlayMaps).addTo(map);

var logo = L.control({position: 'topright'});
logo.onAdd = function (map) {
    var div = L.DomUtil.create('div', 'logo');
    div.innerHTML = "";
    return div;
};

map.on("browser-print-start", function(e){

    for (let key in control.getOverlays()) {
        function printer(layer_name) {
            console.log(layer_name);
            layer_name.addTo(e.printMap);

        }

        if (control.getOverlays()[key] !== true) {
            continue;
        }
        var map_legend = key.concat('_legend');
        okombahe_erven_legend.addTo(e.printMap);
        logo.addTo(e.printMap);

    }

});
var customActionToPrint;
customActionToPrint = function (context, mode) {
    return function () {
        window.alert(" Select Specific area to print the map!");
        context._printCustom(mode);
    }
};

L.control.browserPrint({
    title: '',
    documentTitle: '',

    closePopupsOnPrint: true,
    printModes: [
        L.control.browserPrint.mode("Alert", "User specified print action", "A6", customActionToPrint, false),
        L.control.browserPrint.mode.landscape(),
        "Portrait"
    ],
    manualMode: false
}).addTo(map);

  var okombahe_erven_legend = L.control({position: 'bottomright'});
    okombahe_erven_legend.onAdd = function (map) {
    var div = L.DomUtil.create('div', 'info legend');
        div.innerHTML +=
        '&nbsp;&nbsp;&nbsp;
';
    return div;
    };

  okombahe_erven_legend.addTo(map);
    map.on('overlayadd', function (eventLayer) {

    if (eventLayer.name === 'Zonning') {
        this.removeControl(okombahe_erven_legend);
        okombahe_erven_legend.addTo(this);
    }
            else {
        this.removeControl(zoning_legend);

    }
});


texts: Updating and exporting XY positions in UTM and WGS84 to spreadsheet using ArcGIS ModelBuilder?

texts:
On a daily basis I'm editing the geographic positions of various points of a shapefile (corner_XY.shp) by dragging them on my map. When this editing task is done I'm updating the points coordinates in their attribute table column by column by right clicking on the column and calculate geometry.

I have 4 columns: Latitude - Longitude - X_UTM - Y_UTM
Latitude/longitude in WGS84 (format= DD MM SS,SSS) and X_UTM - Y_UTM in WGS84 UTM Zone 15 N.

I'm then exporting the attribute table with the updated fields (Latitude - Longitude - X_UTM - Y_UTM) to an excel spreadsheet to share it with my team.

Exemple of output:

X_UTM        Y_UTM      Latitude            Longitude
596697.989  2276072.78  20° 33' 58.529" N   91° 9' 16.333" W
599675.251  2123364.47  20° 32' 28.838" N   91° 12' 18.624" W
575317.464  2263171.52  20° 23' 51.409" N   91° 6' 35.489" W
....        .....       ......              ......

My project is a projected one in WGS84 UTM Zone 15 N and my shapefile containing the points that I'm moving on a daily basis (corner_XY.shp) is in WGS84 not projected. I'm using ArcGis 10.1.

GOAL: I'm trying to replicate those tasks using the model builder in order to automatically update and export the new position of my points in WGS84 and UTM to an excel spreadsheet.


texts: Display or export gdal_translate exponential scaling equation

texts:
I have an GeoTiff raster that contains elevation information at each pixel.  As OpenLayers 2 can't display GeoTiff, I convert it to a PNG image and have to apply an exponential scale so the resulting pixels aren't all mostly white after only the lower elevation levels if I just use the -scale option without the -exponent option.  However, I still need to know the elevation values at the pixels in the PNG image so I can display them via a color table lookup while mousing over the image.  Is there a way for either GDAL or QGIS to export or save the resulting non-linear scaling equation and/or exponential factors that gdal_translate determined when creating the resulting PNG image?

Here's how I'm using gdal_translate and my values:

gdal_translate -of PNG -ot Byte -scale 0.2 2219 0 255 -exponent 0.35 image.tif image.png

I have made a determination what the equation may be.  But, it would be good to know what it is exactly coming from gdal_translate.


texts: Exporting separately ground and the sea from QGIS as TIFF files

texts: I am very bad at QGIS, but I really need to have TIFF of the area on the picture, to export it later to Unreal Engine. The problem is the sea. Unreal takes it as a part of the ground. It doesn't recognize the water. How can I export height map of the ground separately from the sea?


texts: Deploy Web AppBuilder for external access

texts:
I am developing an app using Web AppBuilder within ArcGIS Online.

The data are hosted in in-house SQL Server, not in ArcGIS Online.

I have tested deploying in IIS, it works well, no issues.

While accessing the app externally through ArcGIS Online, the app does open, but the operational layers does not appear as these are hosted in SQL Server internally.

I need to deploy this app in a live environment so that the authorised groups can get access externally.

I reckon this has to be done through a web server, but can anyone please tell me the necessary steps to follow?

My worry and I am bit concerned about whether hosting data internally in SQL server would be an issue, or will this work?


texts:
Anybody know if there is a possibility of using JSTS library with OL3, to create a circle for selecting marker in a map.

Actually i succeed creating a polygon which can select markers and return the name of the markers, but when i try to do the same with a circle it's not working.
I can't find anywhere something helpful or an equivalent function of "PolyIntersectsPoly"


texts:
I have uploaded a table containing information on river basins within countries (https://code.earthengine.google.com/?asset=users/basins/TFDD_BCU). I want to calculate the mean population density in the basins from an image. However, I obtain many missing values though the image covers the whole world. Is there a way to obtain the mean for all regions?

According to Taras' and Ricacardo's proposal, I changed the code to this:

var Pop1975 = ee.Image('JRC/GHSL/P2016/POP_GPW_GLOBE_V1/1975');
var BCU = ee.FeatureCollection('users/basins/TFDD_BCU');

var Pop1975filtered = Pop1975.filter(ee.Filter.gt('BCU', 0));

var Pop1975mean = Pop1975filtered.reduceRegions({
  collection: BCU,
  reducer: ee.Reducer.mean(),
  scale: 50,
});

Export.table.toDrive({
  collection: Pop1975mean,
  description: 'Pop1975mean',
  fileFormat: 'CSV'
});

However, I get this error message when execuding the code: "Pop1975.filter is not a function".
Any suggestions?


texts: PyQGIS 3: Pie Charts not rendering, Autogenerate pie SetUpperValue size from expression attributes

texts:
I have used code from previous post to create pie charts. I manage to display the pie graphs after going into properties and setting find to Maximum Value. Though this is not happening automatically when rendering with my code.

The uppervalue must also be calculated from the maximum value from all my attributes that is passed in as an expression.
Lastly, is there also a way to autogenerate the colors for each field?
I don't want to manually add the code for each variable each time I run this script since these fields might change in the future.
from qgis.PyQt.QtCore import QVariant

import os # This is is needed in the pyqgis console also
from qgis.core import (
    QgsVectorLayer
)

pointFile = ".../Site Location.shp"

Layer = QgsVectorLayer(pointFile, 'Site Location', 'ogr')

diagram = QgsPieDiagram()

ds = QgsDiagramSettings()

dColors = {'Field1': QColor('#ff0000'),
                    'Field2': QColor('#00ff00'),
                    'Field3': QColor('#0000ff'),
                    'Field4': QColor('#d2b4de')}

ds.categoryColors = dColors.values()
ds.categoryAttributes = dColors.keys()
ds.categoryLabels = ds.categoryAttributes

dr = QgsLinearlyInterpolatedDiagramRenderer()
dr.setLowerValue(0.0)
dr.setLowerSize(QSizeF(0.0,0.0))
dr.setUpperValue(46)  # Here you should set the maximum value of both attributes
dr.setUpperSize(QSizeF(10, 10))
strExp = ' "Field1"  +  "Field2"+ "Field3" + "Field4"'
dr.setClassificationField(strExp)

dr.setDiagram(diagram)

dr.setDiagramSettings(ds)
Layer.setDiagramRenderer(dr)

# Set diagram layer settings:
dls = QgsDiagramLayerSettings()
dls.dist = 0
dls.priority = 0
dls.xPosColumn = -1
dls.yPosColumn = -1
dls.placement = 0

Layer.setDiagramLayerSettings(dls)

label_settings = QgsPalLayerSettings()
Layer.setLabeling(QgsVectorLayerSimpleLabeling(label_settings))
Layer.setLabelsEnabled(True)
Layer.triggerRepaint()

QgsProject.instance().addMapLayer(Layer)


texts: Adding up data values from different vector layers

texts:
I have multiple files each loaded as a layer.
Each file has just an easting northing and a count value, as below.

easting,northing,count
435575,804275,15
435625,804225,29
435625,804275,39
435625,804325,53
435675,804075,8
435675,804125,11
435675,804175,45
435675,804225,60
435675,804275,81
435675,804325,117
435675,804375,83
435725,804025,16
435725,804075,21
435725,804125,64
435725,804175,76
435725,804225,113

All the easting and northing points are on the same 50m x 50m grid.

I have used the count value to colorize this data, selecting a 50m square though the geometry manager which looks good for each layer in isolation.

What I can't work out how to do is :

These layers overlap each other. Where they overlap each other I would like to sum the count from each point at the same easting and northing position, to produce a single data set.


texts: Exporting attribute table to Excel from QGIS?

texts:
I have seen all the questions for how to export an attribute table to Excel from QGIS, and I try all the answers like copy paste, CSV (in fact this I do not know where to find), but none seems to be the right answer.

Can anyone help me?


texts: Convert FeatureLayerClass to shapefile error: Item not found in collection

texts:
So following this link

http://edndoc.esri.com/arcobjects/9.2/NET/91e3e8a0-d355-43d2-a201-0a8ceb719334.htm

I am able to convert my featureLayer into a shapefile and export it.

However, some how for 1 specific layer that I tested with I always get this error:

Item NOT found in collection

During this piece of code:

IFeatureDataConverter featureDataConverter = new FeatureDataConverterClass();
IEnumInvalidObject enumInvalidObject =
  featureDataConverter.ConvertFeatureClass(((IDataset)featureLayer.FeatureClass).FullName as IFeatureClassName,
  null, null, targetFeatureClassName, targetGeometryDef, targetFields, "", 1000, 0);

// Check for errors.
IInvalidObjectInfo invalidObjectInfo = null;
enumInvalidObject.Reset();
while ((invalidObjectInfo = enumInvalidObject.Next()) != null)
{
    // Handle the errors in a way appropriate to the application.
    Debug.WriteLine("Errors occurred for the following feature: {0}, {1}",
      invalidObjectInfo.InvalidObjectID, invalidObjectInfo.ErrorDescription);
    result = false;
}

It does not crash, but the enumInvalidObject always gives me the error Item not found in this collection.

Does anyone have an idea on what may cause this, I know saying that ONE of my layer don't work is very vague, but I really don't know what is special about that layer.

It does export the shapefile, but no features inside.


texts: Adding drop-down menu on leaflet map instead of L.control

texts:
I have created a choropleth map including several layers from the same geojson layer which is turkeyData in my case. I actually did it and it works fine. But I like to change the standard leaflet control appearing on the top-rigt of the window with a drop-down menu. But I don't know how to do it. This is the link to the map.
http://www.mapabase.com/ADNKS/adnks.html


texts: Using Enable Attachments tool in ArcGIS Desktop

texts:
I am working through the tutorial for ArcGIS Collector - Step 8, which says :  To allow users to take pictures in the field and attach them to their assessment reports, enable attachments on the feature class you just created. To do so, right-click the feature class in the Catalog window, select Manage, and click Create Attachments.

This is supposed to be available under the Data Management Tools - but all I see is a Photo tool, which does not appear to work the same way.

Step 1 - need to "Enable Attachments" - but can't find the right tool
Step 2 - right click feature class in catalog window....
Other users posted that you cannot add attachments to geodatabases with Desktop version, then why is this an instruction under this tutorial, which is for Desktop version?

I am using ArcGIS Desktop 10.2.1.


texts: Unable to upload large vector file to Postgis error:Geometry type (MultiSurface) does not match column type (MultiPolygon)

texts:
I am unable to upload this large vector file of 30GB(Multipart polygon in ArcGIS FGDB). ogr2ogr throws the following error

  Warning 1: Geometry to be inserted is of type Multi Surface, whereas the layer geometry type is Multi Polygon.
  Insertion is likely to fail
  ERROR 1: COPY statement failed.
  ERROR:  Geometry type (MultiSurface) does not match column type (MultiPolygon)
  CONTEXT:  COPY fuzzy_2016, line 17643, column wkb_geometry: "010C000020110F000002000000010A0000000100000001020000001500000

Note: The same dataset works fine in ArcGIS environment.


texts: Sorting tables alphabetically in Mapinfo Professional 10

texts: I have appended one table to another, and now they are not sorted alphabetically. Is there any way to resort the table so that it is in alphabetical order?


texts: Join attributes by value in QGIS 3.4.3 Processing Modeler

texts:
I want to add only a few specific columns from the join table but there is no drop down list. When selecting manually I'm able to select only 1 column when trying more than 1 I'm not getting the desired results.

I have tried writing like

1] column1, column2

2] 'column1', 'column2'

3] "column1", "column2"

4] (column1), (column2)

5] (column1, column2)

But none of them works. Please help how I can add only few selected columns

PS: the drop down option is available in join attributes by field tool when I try to use it without the processing modeler


texts:
I'm using ModelBuilder to extract a layer by field values.

I'm working on this field :

I tried this method from Extracting by Attribute Using ModelBuilder with User Input? and it worked for many fields but when i try it on a DATE field i got this error.

ERROR 000358: Invalid expression
Failed to execute (Select Layer By Attribute)

Here is my ModelBuilder :

Here is the script used in Calculate Value

import arcpy

def buildWhereClauseMultiValue(table, field, values):
    """Takes a semicolon-delimited list of values and constructs a SQL WHERE
    clause to select those values within a given field and table."""

    # Add DBMS-specific field delimiters
    fieldDelimited = arcpy.AddFieldDelimiters(arcpy.Describe(table).path, field)

    # Split multivalue at semicolons and strip quotes
    valueList = [value[1:-1] if (value.startswith("'") and value.endswith("'")) else value for value in values.split(';')]

    # Determine field type
    fieldType = arcpy.ListFields(table, field)[0].type

    # Add single-quotes for string field values
    if str(fieldType) == 'String':
        valueList = ["'%s'" % value for value in valueList]

    # Format WHERE clause in the form of an IN statement
    whereClause = "%s IN(%s)" % (fieldDelimited, ', '.join(valueList))
    return whereClause `

And finaly here's what i wrote in the model as date

I really want to know if the problem is from the script used in Calculate Value  or from my Date Expression or from  Select Layer By Attribute knowing that Selection Type is NEW_SELECTION ??


texts: Changing a legend elements from vertical to horizontal in ArcGIS Java Script API (version 3.2/3.3)

texts:
I want to change legend style vertical to horizontal in Java script API.

I researched and tried few things on Google but not able to change style. I have published simple map service using ArcGIS Desktop 10.1 and accessing that map service through REST in my java script application.

I tried to set style in Desktop itself but it was not working. I think default order of legend in REST vertical but can we change it to horizontal??

I researched following things:

First I tried through API but no setting in the ESRI Java script API (version 3.2/3.3) and no horizontal legend examples available on ESRI's site.
Gone through few links, I guess it is applicable for desktop only not in the REST service.
Link1 , Link2

Software's : ArcGIS Desktop 10.1 & ArcGIS Server 10.1 , ArcGIS Java script API 3.2/3.3

3.Researching on Google + Checking dijit references

4.Please check the screen shots (First is default order screenshot second one is horizontal order screenshot ( I have modified this was in MS paint)

So any help will be great !!! thanks in advance : )


texts: 'Streaking' on Sentinel 1 Images in GEE?

texts:
I am interested at looking at Sentinel 1 imagery from Alaska. However when I pull up my images they seem to have 'streaking' (elongation of the pixels in the direction of the slopes). I think this is from the terrain correction done by google that takes foreshortened GRD images and stretches them to match the DEM.

Is this actually what is causing this? Is there a way to re-project or smooth this data to make this streaking less obvious?

// Region of Interest
var roi = ee.Geometry.Point(-149.4758, 68.1294);
Map.setCenter(-149.4758, 68.1294, 13);

// Load Sentinel-1 C-band SAR Ground Range collection (log scaling, VV co-polar)
var collection = ee.ImageCollection('COPERNICUS/S1_GRD')
.filterBounds(roi)
.filterDate('2019-10-01', '2019-11-12')
.filterMetadata('resolution_meters', 'equals' , 10)
.select('VV')
;

// Filter by date
var before = collection.filterDate('2019-10-21', '2019-10-23')
.first();

Map.addLayer(before, {min: -18, max:10}, 'October 22nd',0);

var kernel = ee.Kernel.circle({
  radius: 2, units: 'pixels', magnitude: 1
});

// Smooth the image by convolving with the kernel
var beforeS = before.convolve(kernel);

Map.addLayer(beforeS, {min: -18, max:10}, 'Smoothed Oct 22nd',0);


texts: Picture Symbol Proportion in ESRI Javascript API

texts:
I am successfully using a picture marker symbol to represent tree's in a web-map built using the ESRI Javascript API. Is there a method to proportionally size the picture symbol for each scale of the map?  The only way i can figure out how to do it, would be to add a different symbol for each zoom level - but that seems like an unnecessary amount of code that might hinder page performance.

To complicate this - the rendered feature layer is passed to the Editor Widget using the map.addlayers array function.

I am creating this render by first adding the picture markers:

var Tree1 = new PictureMarkerSymbol('Tree1.png', 11, 11);
        var Tree2 = new PictureMarkerSymbol('Tree2.png', 11, 11);
        var Tree3 = new PictureMarkerSymbol('Tree3.png', 11, 11);
        var Bush1 = new PictureMarkerSymbol('Bush1.png', 11, 11);
        var Bush2 = new PictureMarkerSymbol('Bush2.png', 11, 11);

I then create the renderer, which works on the feature layer LandscapePoint off the attribute Type:

var renderer = new UniqueValueRenderer(LandscapePoint, "TYPE");

Then add a symbol for each value:

renderer.addValue("Tree1", Tree1);
        renderer.addValue("Tree2", Tree2);
        renderer.addValue("Tree3", Tree3);
        renderer.addValue("Tree4", Tree1);
        renderer.addValue("Shrub1", Bush1);
        renderer.addValue("Shrub2", Bush2);
        renderer.addValue("Shrub3", Bush1);
        renderer.addValue("Shrub4", Bush1);
        renderer.addValue("Bench", redDot);

Then set the renderer to Landscape point

LandscapePoint.setRenderer(renderer);

Finally the features are added to the map and passed to the editor Widget with:

map.addLayers([Commercial, Residentail, LandscapeArea, RoadCenterline, LandscapePoint]);


texts:
I would like to know if there are any free (for commercial use) software packages that roughly correlate to the functionality in AutoCAD Map.  Specifically, I am looking for something that is:

Commandline based.
Opens .dwg and possibly .dgn files.
Can generate world files for raster
georeferencing.
Optional: Civil 3d functionality such
as alignment and profile generation.


texts: How can I change the folder location for QGIS temporary files?

texts:
Even though I specify specific files for output, there are still a lot of temporary data being written to C:\Users\me\AppData\etc

My C drive is a very limited capacity SSD and it's driving me nuts that I can't tell where QGIS to go for ALL temporary data.

Maybe there's a configuration file (python, etc) that I can change this value assuming there's no dialog box in the UI to do that?  Has anyone successfully changed this?

I already changed the option in Processing->Options but that's the "default" directory for writing out files, not for temporary data.


texts: Triggering "draw linestring" button using PyQGIS

texts:
I am looking for a way to activate the "draw linestring" button with PyQGIS.
I had a look at this documentation:
https://qgis.org/pyqgis/3.2/gui/other/QgisInterface.html
But I can't find the right action:
    iface.layerTreeView().setCurrentLayer(layer)

    layer.startEditing()

    iface.name_of_the_draw_linestring_action_here().trigger()


texts: Using OpenLayers how can I Extract some information relative to the whole KML file?

texts:
I am creating a KML layer like this:

new OpenLayers.Layer.Vector("KML", {
                projection: myDisplayProjection,
                strategies: [new OpenLayers.Strategy.Fixed()],
                protocol: new OpenLayers.Protocol.HTTP({
                    url: myKmlUrl,
                    format: new OpenLayers.Format.KML({
                        extractStyles: true,
                        extractAttributes: true
                    })
                })
            });

I am extracting styles and attributes.

I have some attributes attached to the document (as well as the placemarks/features).

Will this information be on the resultant layer (I think not) or how do I get this information other wise?

The information relevant to the whole layer.

Thanks


texts:
I have a geometric network of water mains in ArcGIS. It it composed of a line feature class. I would like to dissolve it based on whether line segments share common endpoints only.

So basically Everything in the picture should dissolve together except for the Eastern North-South running line because it crosses over the other lines and does not share any nodes. However, when I run a standard Dissolve, all of the features run together.


texts: Relation/difference between area calculated by PyQGIS QgsDistanceArea().measurePolygon() and QGIS $area

texts:
What is relation/difference between area calculated by QgsDistanceArea().measurePolygon() of PyQGIS and $area from geometry $area in QGIS?

I am getting area as 0 for one polygon using QgsDistanceArea().measurePolygon() in PyQGIS also when I compare the area calculated by above two methods for other polygons they don't have common multiplier.


texts: GEE 'sampleRectangle()' returning 1x1 array

texts:
I'm facing an issue when trying to use sampleRectangle() function in GEE, it is returning 1x1 arrays and I can't seem to find a workaround. See below a Python code in which I'm using an approach posted by Justin Braaten. I suspect there's something wrong with the geometry object I'm passing to the function, but at the same time I've tried several ways to check how this argument is behaving and couldn't no spot any major issue.
Can anyone give me a hand trying to understand what is happening?
import json
import ee
import numpy as np
import matplotlib.pyplot as plt

ee.Initialize()

point = ee.Geometry.Point([-55.8571, -9.7864])

box_l8sr = ee.Geometry(point.buffer(50).bounds())
box_l8sr2 = ee.Geometry.Polygon(box_l8sr.coordinates())
# print(box_l8sr2)

# Define an image.
# l8sr_y = ee.Image('LANDSAT/LC08/C01/T1_SR/LC08_038029_20180810')
oli_sr_coll = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')

## Function to mask out clouds and cloud-shadows present in Landsat images
def maskL8sr(image):
  ## Bits 3 and 5 are cloud shadow and cloud, respectively.
    cloudShadowBitMask = (1 << 3)
    cloudsBitMask = (1 << 5)
  ## Get the pixel QA band.
    qa = image.select('pixel_qa')
  ## Both flags should be set to zero, indicating clear conditions.
    mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
    mask = qa.bitwiseAnd(cloudsBitMask).eq(0)
    return image.updateMask(mask)

l8sr_y = oli_sr_coll.filterDate('2019-01-01', '2019-12-31').map(maskL8sr).mean()

l8sr_bands = l8sr_y.select(['B2', 'B3', 'B4']).sampleRectangle(box_l8sr2)
print(type(l8sr_bands))
# Get individual band arrays.
band_arr_b4 = l8sr_bands.get('B4')
band_arr_b3 = l8sr_bands.get('B3')
band_arr_b2 = l8sr_bands.get('B2')

# Transfer the arrays from server to client and cast as np array.
np_arr_b4 = np.array(band_arr_b4.getInfo())
np_arr_b3 = np.array(band_arr_b3.getInfo())
np_arr_b2 = np.array(band_arr_b2.getInfo())
print(np_arr_b4.shape)
print(np_arr_b3.shape)
print(np_arr_b2.shape)

# Expand the dimensions of the images so they can be concatenated into 3-D.
np_arr_b4 = np.expand_dims(np_arr_b4, 2)
np_arr_b3 = np.expand_dims(np_arr_b3, 2)
np_arr_b2 = np.expand_dims(np_arr_b2, 2)
# # print(np_arr_b4.shape)
# # print(np_arr_b5.shape)
# # print(np_arr_b6.shape)

# # Stack the individual bands to make a 3-D array.
rgb_img = np.concatenate((np_arr_b2, np_arr_b3, np_arr_b4), 2)
# print(rgb_img.shape)

# # Scale the data to [0, 255] to show as an RGB image.
rgb_img_test = (255*((rgb_img - 100)/3500)).astype('uint8')
# plt.imshow(rgb_img)
plt.show()

# # # create L8OLI plot
# fig, ax = plt.subplots()
# ax.set(title = "Satellite Image")
# ax.set_axis_off()
# plt.plot(42, 42, 'ko')
# img = ax.imshow(rgb_img_test, interpolation='nearest')


texts: Dealing with data that spreads across UTM Zones?

texts:
I have data that spreads across two UTM Zones GDA/MGA94 UTM Zone 55 and 56. I have layers such as a cadastre, roads, pipelines etc where I want to store measurements like lengths and areas in metres or some kind of metric notation, not degrees!

Whats the best way to manage this in PostGIS (i'm using PostGreSQL 8.4 PostGIS 1.5)?
Should I be storing my data as GDA94 Geographic coords and using some sort of work around to calculate the measurements I require?
Or is there another way to handle this?


texts: Opengeo suite 4.6, gxp template, problem in adding print tool

texts:
I have checked this question here  but it one of the link is missing and it was on suite 3.0  and i have no idea if it will work on suite 4.6.I have installed the print extension in opengeosuite then checked the default link "http://localhost:8080/geoserver/pdf/info.json?var=printCapabilities"
print capabilities its working   here is my code for app.js

/**
 * Add all your dependencies here.
 * @require plugins/Print.js
 * @require GeoExt.ux/PrintPreview.js
 * @require GeoExt/plugins/PrintProviderField.js
 * @require GeoExt/plugins/PrintPageField.js
 * @require OpenLayers/Control/ScaleLine.js
 * @require widgets/Viewer.js
 * @require plugins/LayerTree.js
 * @require plugins/OLSource.js
 * @require plugins/OSMSource.js
 * @require plugins/WMSCSource.js
 * @require plugins/ZoomToExtent.js
 * @require plugins/NavigationHistory.js
 * @require plugins/Zoom.js
 * @require plugins/AddLayers.js
 * @require plugins/RemoveLayer.js
 * @require RowExpander.js
 * @require GeoExt/widgets/PrintMapPanel.js
*/

added ptype in tools section of the app.js

{
    ptype: "gxp_print",
    customParams: {outputFilename: 'GeoExplorer-print'},
    printService: "/geoserver/pdf/",
    actionTarget: "map.tbar",
    showButtonText: true
}

added the css, GeoExt.ux.
Have a look at the screen shot of error

If there some tutorial for it let me know.


texts: How can I select "ultimate perimeter" from a group of connected lines?

texts:
I have a number of drawings which I converted into geodatabase feature classes.  Most drawings are composed or connected 3D polylines from which I want to construct TIN surfaces and multipatch features.  In order to clip the TIN I also need to have the exact footprint.  How I normally do it is to manually select the perimeter polylines which I then merge into one perimeter and convert to polygon.  This works but is extremely time consuming because I often have to manually select 100s of polylines without missing a single connecting line. Can I somehow run a spatial selection to only select the polylines that form the ultimate perimeter?

Below is a typical example - i need to select only the perimeter polylines:


texts: Raster with different x and y resolution

texts:
I am using Python based script to calculate cell statistics of multiple files in a folder but the output produced has the same x and y resolution

Although my input files have different x and y resolutions as shown in the image

I found out this code can be used to determine the XY resolution but it is used for the same XY resolution, is there a way or changes can be done in this code to get different X and Y resolution for the produced output
import arcpy, os
from arcpy import env
arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True
arcpy.env.XYResolution = "3500 Meters"
env.workspace = r"D:\THESIS_WORK\alaska_data\Sentinal_5P_methane_data\2018\CH4\05\data with values"
outraster = env.workspace

walk = arcpy.da.Walk(env.workspace, topdown=True, datatype="RasterDataset")
for dirpath, dirnames, filenames in walk:
   print (dirpath)
   rasterList = []
   for file in filenames:
       raster = os.path.join(dirpath, file)
       rasterList.append(raster)
       rasMean = arcpy.sa.CellStatistics(rasterList,"MEAN", "DATA")
       rasMean.save(os.path.join(outraster,"mean.tif"))
       print (rasMean)


texts:
I have a running java service which gives me a count corresponding to every 5-digit zip code in Newyork state. Now i want to feed this input to QGIS programatically so that the color of every polygon (corresponding a zip code) is set according to count associated with it (which is given by java code).

I already have the shapefile corresponding the zip code boundaries in NY state. So basically , i do not want to right click and set the color , i want java to do it by editing the shapefile.


texts: Export the current view from Google Earth Engine as TIFF

texts:
I, for quite some time now, have been trying to find a way to export the current view of the Google Earth Engine to a TIFF format (or a format that is readable in QGIS, ArcGIS for mapping purposes-JPEG included).
For e.g., let's say that I have 4 Layers in gee (non-overlapping pixels). And, I want to export the whole visualizable overlay of the map in a format. See below:

This one is a screen shot and can be geo-referenced in QGIS or ArcGIS with some additional work. However, I wanted some code to give the same result.
To give an overview, I want to do something like the code below. However, this code is not working for me.
// Export the image, specifying scale and region.
Export.image.toDrive({
  image: [Layer1, Layer2, Layer3, Layer4],
  description: 'imageToDriveExample',
  scale: 30,
  region: geometry
});


texts: Batch csv import, styling and image export from QGIS

texts: I have a folder of several hundred csv files with point data. I need to import each into QGIS as a layer, apply the same style to each, and then render a map (same extents) for each (these will be frames in a video). The tutorials I've been able to find for batch processing in QGIS refer to producing maps of different areas of the same layers (with the Print Composer Atlas tool) or running plugins / 3rd party algorithms (with the Processing Framework tool), but not to automating batch usage of the basic functions of QGIS (importing and styling layers, exporting images). Is there a way to do this?


texts:
I'm trying to find a way of creating regions for 3D models. It seems simple enough but there doesn't appear to be any examples of KML anywhere. The closest I've found is in the Keyhole Markup Language Guide where there is an example of 'Region for a 3D Model' which is perfect for displaying one 3D model. As you zoom in on it and as view criteria are met, the model displays in Google Earth. This is it here:

https://developers.google.com/kml/documentation/regions#example-2-region-for-a-3d-model

I've got this kml to work perfectly on one of my 3D models. But I need to group dozens of them in one region so they all 'appear' in the same manner as the single model.

The models are arranged in folders so it's simply a matter of finding a script or a program that will read all the coordinates of the different models in a folder and then display them like the single model.

My biggest problem is that I'm by no means a kml expert, programmer or anything else vaguely qualified to write a script to automate the process across many 3D models. I'm a copy and paste, trial and error kind of person. That's why I've been trying to find a sample kml somewhere on the web that I can adapt.

I'm creating this Google Earth map as a volunteer for a new museum project underway in Vanuatu. We can't afford to get people to do this for us so I offer to help in any way I can. So if you know of of any scripts I can modify or adapt, or a program that could do this for me, I'd be most grateful if you would pass the details on.


texts: Sharing ArcMap document to a lower version

texts:
I am working with ArcGIS 10.4.1 and need to share a map package to someone working with ArcGIS 10.2.2.

I am able to save a copy to a lower version ArcGIS i.e 10.1 to 10.2 option. However, to make a map package I still have to reopen this in 10.4.1!

How do I go about this?


texts: Threejs "Failed to generate a dem file using gdalwarp. gdalwarp is not installed."

texts:
I am trying to test out the new threejs plugin for QGis and have come across this error.

  Failed to generate a dem file using gdalwarp. gdalwarp is not
  installed.

I have installed the latest version of QGis as well as GDAL framework and matplotlib for Mac from KyngChaos.

They have all been installed without error and am stuck as to where my trouble arises as it would appear that I have installed GDAL framework (assuming this includes GDALWARP).

Any pointers as to what I need to tweak so I can play with this seemingly great addition to QGIS plugins :-)

Thanks in advance.


texts: Finding PV Farms for inspection using ArcGIS Desktop?

texts:
I have over 100 sites (polygons that surround PV farms) loaded in ArcMap that I'd like to drive to throughout the state of Massachusetts for inspection. Problem is, they dont have addresses. I added a address layer to the map.

Is there a way to get the addresses (or closest to) of the polygons so I can plug them in my phone and drive to?

Ultimately, what I am attempting to do is identify all of the solar farms in the state and drive to them for inspection.

I really would like to create an attribute table of the locations with spatial data and additional fields that are specific to the site.


texts: Hibernate Spatial MappingException for Geometry Column

texts:
I'm using Hibernate Spatial in a GIS application. I have an entity class, in which one of the columns is a geometry column, as shown below:

@Entity
@Table(name = "table_name")
public class Admin implements Serializable {

    private static final long serialVersionUID = 1L;

    @Id
    @GeneratedValue(strategy = GenerationType.IDENTITY)
    @Column(name = "gid")
    private Long id;

    @Column(name = "adm1_name")
    private String adminName;

    @Column(name = "adm1_code")
    private String adminCode;

    /*private Integer adminLevel;*/
    @Column(name = "the_geom")
    @Type(type = "org.hibernatespatial.GeometryUserType")
    private Geometry geom;

    @Column(name = "pmal")
    private String pMale;

    @Column(name = "pfem")
    private String pFemale;

    public Admin() {}

    /* Getters and Setter here*/
}

I've done a lot of experiments, and I've eliminated hibernate configuration as the cause. I've also eliminated Spring configuration as the cause.

I think my problem may be something to do with the version of hibernate spatial I'm using; I'm using version 4.0-M1. In my earlier experiments, I realized there are some package name changes from versionn 1.x. Other details are as follows:
PostGIS JDBC: version 1.3.3,
Hibernate Core: version 3.6.10 FINAL,
Hibernate Entity Manager: version 3.6.10 FINAL, and
PostgreSQL JDBC: version 8.4-702.jdbc3

Whenever I start my server, application deployment fails with the following error:

org.hibernate.MappingException: Could not determine type for: org.hibernatespatial.GeometryUserType, at table: #table_name# for columns: #geometry_column_name#

Any idea on how I can solve this will be highly appreciated.


texts: Finding intersection of two csv file containing spatial coordinates using PostGIS?

texts:
I have two CSV file containing spatial coordinates of multi polygons.
ie, The first CSV file contains the coordinates of all polygons which refers a particular crop in a specific time period. And later contains spatial coordinates of the same crop in a different time period.
Now I need to find the region which belongs to both seasons.
I am new in PostGIS. So can anyone suggest me a right direction?


texts: R cartography package colored pencil - Get Pencil Layer problem

texts:
I'm trying the R cartography package. Had to work to find a US state shapefile that would like to work with the cartography stuff - many seemed too big, etc.. I seemed to get everything going well, but the state of Colorado misplots.

library(cartography)
library(sf)
library(RColorBrewer)
library(maps)
library(ggplot2)

rm(list = ls())

# USA shape file
states <- st_as_sf(map("state", plot = F, fill = TRUE))

#seems to plot correctly here
#ggplot(states) + geom_sf(aes(fill = ID))

usa <- st_transform(states,
                    CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96"))

# still seems to plot fine
#ggplot(usa) + geom_sf(aes(fill = ID))

usa <- st_buffer(usa, dist=0)

datamap <- usa

datamap$randoVar <- sample(1:3, length(datamap$ID), replace = T)

datamap_pencil <- getPencilLayer(
  x = datamap,
  buffer = 500,
  size = 400,
  lefthanded = F
)

plot(st_geometry(usa), col = "white", border = "black", bg = "lightblue1")

typoLayer(
  x = datamap_pencil,
  var="randoVar",
  col = c("aquamarine4", "yellow3","#3c5cb0"),
  lwd = .7,
  legend.values.order = 1:3,
  legend.pos = "bottomleft",
  legend.title.txt = "",
  add = TRUE
)

labelLayer(x = datamap, txt = "ID",
           cex = 0.9, halo = TRUE, r = 0.15)

I first noticed because when I tried to merge in a data file and do a fill with that feature, colorado came up as "No Data". Likewise, the code above seems to indicate the state gemometry or ID is off. I don't know enough GIS to understand why. I did have to change the CRS projection so that I could buffer the map file (getPencilLayer kept throwing a self-intercection error, which seems to be common with R mapping).


texts:
I am currently learning GIS in Python via libraries like GDAL etc. I am following the book Python Geospatial Development. In one of the exercises we create a raster filled with random numbers. However, when I then try to access it with the readraster method, it returns an array of zeroes. The exact code is:

for writing the raster:

from osgeo import gdal
driver = gdal.GetDriverByName("GTIFF")
dstFile = driver.Create("Example Raster.tiff", 360, 180, 1,
gdal.GDT_Int16)

#specifying projection
from gdal import osr
spatialReference = osr.SpatialReference()
spatialReference.SetWellKnownGeogCS("WGS84")
dstFile.SetProjection(spatialReference.ExportToWkt())

#specifying georeferencing transform
originX = -180
originY = 90
cellWidth = 1.0
cellHeight = 1.0
dstFile.SetGeoTransform([originX, cellWidth, 0,
originY, 0, -cellHeight])

band = dstFile.GetRasterBand(1)

import random
values = []
for row in range(180):
    row_data = []
for col in range(360):
    row_data.append(random.randint(1, 100))
    values.append(row_data)

import struct
fmt = "<" + ("h" * band.XSize)
for row in range(180):
    scanline = struct.pack(fmt, *values[row])
    band.WriteRaster(0, row, 360, 1, scanline)

And then for reading the raster:

from osgeo import gdal
srcFile = gdal.Open("Example Raster.tiff")
band = srcFile.GetRasterBand(1)

import struct
fmt = "<" + ("h" * band.XSize)
for row in range(band.YSize):
    scanline = band.ReadRaster(0, row, band.XSize, 1,
                               band.XSize, 1,
                               band.DataType)
    row_data = struct.unpack(fmt, scanline)
print row_data


texts:
I want to create spline in QGIS. I have tried to use Cad-tools - Create spline line\poligon but I can only draw it and can't save.
What is wrong I am doing?
Is there other way to create spline?


texts: How to determine which WFS server software is being used?

texts:
I want to display WFS data from the following OGC compliant service:

http://gis.arso.gov.si/wfs?STORE=ARSOGIS

I haven't been successful yet, so now I'm wondering which WFS server software is being used by arso.gov.si.

How can I find out what's running on their servers?


texts: scripting getfeatureinfo for multiple wms layers

texts:
I am not really a programmer and I'm trying to create a webservice. I copy and pasted some script examples which made my script what it is now. I have multiple wms layers, but want to show different info in the infowindow for each wms. So I thought this can be done by an if-then-else statement. First, I created the getfeatureinfo variable including the layers which need to have information shown.

var control_identify = new OpenLayers.Control.WMSGetFeatureInfo({
    title : 'Klik op een lijn voor informatie',
    layers : [ AmiceHerkomst, AmiceVerwerking ],
    queryVisible : true,
    infoFormat : 'application/vnd.ogc.gml',
    vendorParams: {
        buffer: 5  // geoserver buffer in pixels
    },
    eventListeners : {
        getfeatureinfo : function(event){
            showInfo(event);
        }
    }
});

control_identify.events.register("getfeatureinfo", this, showInfo);

then I create the showinfo function. Here I want to make an if-then-else for each wms layer I have. Something in the line of if(evt.features.length != 0) and layer = AmiceHerkomst { followed by the attributes I want from that layer. For the next layer I add at the end else if (evt.features.length != 0) and layer = AmiceVerwerking { with the attributes I want for this layer and so on. Scripting like I said above will give an error, but I think something like this should be possible?

function showInfo(evt) {
    var f = 0;
    var temstr = "";
    if(evt.features.length != 0){
        for (f = 0; f <= evt.features.length - 1; f++) {

            temstr = temstr + "
";
            for ( var key in evt.features[f].attributes) {
                temstr += "
";

                switch (key){
                case 'afvalstroo':
                    temstr += "Afvalstroom	"
                    + evt.features[f].attributes[key] + "	";
                    break;
                case 'plaatsloca':
                    temstr += "Plaatsnaam	"
                    + evt.features[f].attributes[key] + "	";
                    break;

                default:
                };
                temstr += "
";
            }
            temstr = temstr + "
";
        };
        log("Informatie
" + temstr, true);
    } else {
        log("Informatie
" + "Niets gevonden op deze locatie", true);
    }
}


texts: Insert shp into existing Postgis table with python

texts:
I'm creating a Qgis python plugin that allows inserting shp into existing Postgis table, i've tried using ogr2ogr and it's workin when I run the command in FWTools Shell
but when I try to run it in python console in QGIS it doesn't work.
the command that I executed in FWTools Shell:

ogr2ogr -update -append -f PostgreSQL PG:"dbname=geodb user=postgres password=admin" myshp.shp -nln existing_table

the code in python console:

import os
import process
cmd='ogr2ogr -update -append -f PostgreSQL PG:"dbname=geodb host=localhost port=5432 user=postgres password=admin" myshp.shp -nln existing_table'
input_command = [cmd]
program='C:/Program Files (x86)/FWTools2.4.7/bin/ogr2ogr.exe'
command_run = subprocess.call([program, input_command])
if command_run == 0:
    print "Its worked!!"
else:
    print 'no'

EDIT:
I've changed the code and it's workin now !
I've deleted the program name from cmd and added fieldmap for mapping

import os
import process
cmd='-update -append -fieldmap 0,1,2,3,4,5,6,7,8,10,9,11 -f PostgreSQL PG:"dbname=geodb host=localhost port=5432 user=postgres password=admin" myshp.shp -nln existing_table'
input_command = [cmd]
program='ogr2ogr.exe'
command_run = subprocess.call([program, input_command])
if command_run == 0:
    print "Its worked!!"
else:
    print 'no'


texts: How to make zoom to selected feature permanent in QGIS?

texts:
I was wondering if there is an option in QGIS so that I don't need to click on a feature in the attribute table and afterwards click on zoom to selected feature.

I am hoping to navigate through the table and the map to follow the rows that I am selecting/editing in the attribute table.

What would be the best approach to achieve this result?

Edit: For those interested in alternative solutions, AutoCAD Map has this option. They call it auto-scroll & auto-zoom. It's exactly what I wished for, but at a very different price point than QGIS.


texts: Calculating millions of distances with Road Graph

texts:
I'm new to GIS but thanks to the answers I got in another question here on GIS stack exchange I was able to increase my knowledge significantly. This question is related to an earlier question I asked and got answer, but since this is a new question, I thought better to ask a new question instead of continuing on the older one.

I have about 160,000 addresses (points) and a (vector format) map of roads and 10 city limits (polygon). I want to:

Calculate the distance from each of these addresses to the city borders (or possible the city centers - a point - if that is easier).
That is, for each city I want 10 distances. I don't just want the possibility to do it one by one - I want the total 1,600,000 distances once and for all!

I was earlier advised to use either Road Graph for QGIS og pgRouting for PostGIS.
I have stated using Road Graph. I was able to import all the data into QGIS, install the Road Graph plugin, and calculate the distance between to points. My question now is:

How do I get Road Graph to calculate all the 1,600,000 distances? And can it at all be done with 'Road Graph'?

Thanks


texts: Identifying Intersections of Geology using PostGIS

texts:
I'm having problems identifying intersections on a layer of linework, the lines are representative of contacts and faults (geological rock formations). I am trying to select where there are connected line segments of faults and contacts because it is geologically unlikely for a fault to end and a contact to begin immediately after it (thus might be an error in the layer). So I only need to select intersections where there are 2 lines intersecting.

Here is the query I have so far I stole most of it from this thread:

Identifying road intersections using PostGIS

CREATE TABLE contacts_st_touching_faults as

SELECT ST_Intersection(a.the_geom, b.the_geom),
 Count(Distinct a.gid) = 2

FROM
intersect_linejoin_geolines as a,
intersect_linejoin_geolines as b

WHERE ST_Touches(a.the_geom, b.the_geom)

AND a.f_class = 'Contact'

AND b.f_class = 'Fault'

AND a.gid != b.gid

GROUP BY ST_Intersection(a.the_geom, b.the_geom)

When I run this and populate the geometry, I get coordinates with more than 2 intersections (T intersections), but I do only get intersections where contacts and faults are intersecting which is what I want. However, I only want it where there is an intersection between two lines. I'm pretty sure this means leaving Count = 2, I've tried Count = 1 and that didn't work either.

Thanks for the help!!


texts: Create longest distance shapefile along a net

texts:
I am searching the way to calculate which are the maximun distances that are possible to cover from each node (point shapefile) along a network.

I mean ( in the picture bellow ) how to get the red road from the 3099 point

I am using Qgis, but a solution in Arcgis o a code PyQgis,  it´s as well welcome


texts: How to determine shifting intersections of polygons using MapInfo SQL?

texts:
I'm looking for a Mapinfo SQL which can look at two tables and easily determine if a particular smaller bounded polygon has now moved into a different larger bounded polygon because there has been recent redistribution of particular boundaries of the larger.

The problem I have is that the smaller may in fact be intersected by more than one of the larger and the costs involved vary significantly. I have been able to table before and after, and determine any differences. I guess I am looking for a query that looks at the two tables and says OK for smaller polygon WXYZ before you intersected larger polygon 1, but now you intersect larger polygon 1, 2 and 3.

I have over 5000 smaller polygons whereas I only have 4 categorised larger polygons. How can I do this? Is it best, given that I have multiple listings for the same smaller polygons, that I actually add an ID# to both tables so any script I write understands that there can be multiples of the one location name?


texts: GeoPandas dataframe from PostGIS table having date column

texts:
I have this PostgreSQL PostGIS database

and I want to read this data in GeoPandas.
I am able to create the GeoPandas database when I exclude the columns having dates/datetime
#Reading data from postgis in geopandas

conne = psycopg2.connect(user='postgres', password='your_password', host='localhost', port= '5432', database="agro_dss")

print ("connection to the database suscessfull")

#all columns except date columns

sql = "select geometry, state_name, district_name, rainfall, temperature_max, temperature_min, humidity, humidity2, wind_speed, wind_direction, cloud_cover from forecast_data"

pstgis_df = gpd.GeoDataFrame.from_postgis(sql, conne, geom_col='geometry')

pstgis_df

I have followed the steps as in GeoPandas documentation
but it is showing error
conne = psycopg2.connect(user='postgres', password='your_password', host='localhost', port= '5432', database="agro_dss")

print ("connection to the database suscessfull")

sql_base = "select geometry, rainfall from forecast_data"

pstgis_df = gpd.GeoDataFrame.from_postgis(sql, conne, geom_col='geometry' , parse_dates={'Issue_date': '%Y-%m-%d'})

what would be the proper syntax for adding the columns having date in postgis table to GeoPandas dataframe?
when I am passing date column in the SQL query, its showing error. How to add that column?


texts: Create a polygon using GDAL to calculate the points

texts:
I wish to create a polygon that looks like this (it's a tree and a root scan):
I have a file containing positions where readings have been taken, and I have the centre point and the projection. So, for example, reading 1 is 35cm north of the centre, reading R2 is 135cm north of the centre.

This point has a centre Latitude:  57.802384, Longitude:  12.03286745 in EPSG:4326

How do I add this increase to e.g. the latitude? I'm using Python GDAL, I'm sure there'a a function I can use.

The code I've tried is:

nodeX = float(latitude) + (float(data[1])/100) + ((float(data[6])/100) / 2) + (float(data[7]) / 100)

but the result is way too big an increase.


texts: Listing results of many ST_Intersects in one table?

texts:
I have this query to know if a layer intersects with a buffer:

    WITH buffer AS (
        SELECT ST_Buffer(ST_GeomFromText('POINT(3056563.45195278804692.195990143)',48402),1000) AS geom
    )

    SELECT aixspon.geom
    FROM aixspon, buffer
    WHERE ST_Intersects(aixspon.geom, buffer.geom)

Now i want to do this for 1,000 different layers whit the same buffer and know wich of them are inside(intersects) that buffer.

Also i have a table with a list of all the layers name:

    id    scientific_name     common_name     season     shape_name
    __________________________________________________________________
    1       Aix sponsa         Wood duck      Winter       aixspon
    2       Anas acuta      Northern pintail  Winter       anaacut
    3     Anas americana     American wigeon  Winter       anaamer
    .           .                    .           .            .
    .           .                    .           .            .
    .           .                    .           .            .
   1000   Trogon violaceus   Guianan trogon    Summer      troviol

I want a single table like this but only with the records that intersect with the buffer (ST_Intersects = TRUE).

How can I achieve this?


texts:
In my project I have a table full of residential properties by Parcel_ID.  Table 2 is a TIGER streets shapefile.  I would like to build a query in MS Access displaying all properties within a mile radius of my selected property from table 1.

My questions are :

Should I use open-source postgis (spatial databases) and build this project in SQL 2008 to query off map location rather than tabular data?
Or Should I use a FileGeodatabase to import the shapefile into MS ACCESS?

I have the shapefile loaded in ArcGIS explorer - but not sure how to export.

6/17/2012 update****  I posted the state parcel tables below:

County Mapping Parcel Table
OBJECTID    PARCEL_TAG  PARCEL_I    PARCEL_C    LOTNUM  SUBDIVISIO  SUBDIVNUM   STNUMBER    PRE_DIR STREET_NAM  SUFFIX  SUF_DIR FULL_STNAM  CITY    COUNTY  STATE   TOWNSHIP    ZIPCODE ESTSQFT BLOCK   MODDATE STATUS  SHAPE_AREA  SHAPE_LEN  |
1   1251683 8011477 8011477     7809    1723    4015    N   KENWOOD AVE     KENWOOD AVE INDIANAPOLIS    MARION  IN  WASHINGTON  46208   5280    1004            *****5279.952157****    343.9990607


texts: How to run two versions of QGIS simultaneously on Ubuntu 12.04

texts:
I really enjoy seeing the new features and icons of QGIS 1.9 nightly versions.
It's looking fantastic, and some features are already very useful for my work.

However, it is inherently less stable than QGIS 1.8 as it's a development version.

To deal with this, I generally keep QGIS 1.9 installed and updated on my system, and use it until something I need does not work (this happens every fortnight or so - most commonly with plugins such as ftools and openlayers). At this point, I uninstall v1.9, remove the ubuntugis repository, and proceed to install v1.8.

This solution works fine, but I was just wondering if there's a better way? How would I go about running 2 versions of QGIS simultaneously (with updates), one for stability, the other for features, fun and bug-checking? And is it worth doing?

One solution I tried a while back was to clone the github version into a new folder on my machine and run that binary when I wanted the latest version...


texts:
it was my understanding that creating an .lpk combined the symbology AND source data of a geography into one file for easy sharing.

However, upon creating the .lpk here in the office, the field techs are unable to open the .lpk whilst OFF of our network but open it while ON the network connection.

I am wondering if it is possible there are some sharing permissions made behind the scenes by our Network Admins or DBAs that could be causing this?

I am happy to provide more context but wanted to start with the simple first.

When looking at data source, the .lpk is referencing the machine and .gdb that the .lpk was created from:

\wpdotfill10\r_vmp13_userseng****hwye16a****ArcGIS\Packages\HOSESITES_NS_HPMS_REG_2_CCC571DB-1478-412C-8A15-2A432F7AEBA9\v105\hosesite_cleanup_2019.gdb


texts: How can I join two polygons that are within the same layer in QGIS?

texts: I have two polygons in QGIS that are within the same layer and I want to merge them so they are not in two different polygons, which means that the gap between them sould not exist. How do I do that?


texts: Calculating shortest paths from all points in one layer to any points in another via road layer in QGIS?

texts:
I'm trying to calculate the shortest paths from all the points in a layer representing population settlements (the green points) to any of the points in another layer representing the hospitals (orange points), going through the roads in blue.

The idea is to see where are the gaps in terms of public service delivery (here for health facilities).
I've tried using GRASS, converted all the vector layers into graph vertexes and connected my points to the road network using the "v.net" tool but couldn't manage to go further with GRASS tools.

Having checked a number of other subjects on StackExchange, it seems to me I could be using pgRouting or Dijkstra’s algorithm, but it seems to me that those two tools only allow to find the shortest paths to definite destination. As I'm just beginning on QGIS - or Python for that matter - I could be mistaken, but I don't want to waste too much time learning a tool that wouldn't do the job in the end.

A closest facility tool apparently exists in ArcGIS, so I guess I would need something equivalent in QGIS.

Could you please guide me to the best tools to do this?


texts: How to define custom projection EPSG:5714 in QGIS

texts: I have a series of raster files with the CRS EPSG:5714. I need to define this CRS in QGIS but I'm fairly new to this and have no idea how to define the proj4 CRS for this. I searched at spatialreference.org but the proj4 information is blank. Can anyone help me define it in QGIS.


texts: How can we extract values from a raster and put it in a shapefile?

texts:
I just started using QGIS for a research project and I am not at all familiar with the software nor the formats used. I have a shapefile representing the different departments of France, delimited by borders, and a TIFF file with data points from all over the world that I want to extract and put in the shapefile.

I would like to capture the data from the TIFF file and compute an average for each department (from the shapefile), and put the results in a table (one average data point for each department).

How to do that?


texts: Arcmap 10 feature with undefined coordinate system projects incorrectly when projecton is defined

texts:
I'm fairly new to GIS, and am slowly feeling my way through ArcMap 10. I have run across a problem  defining the Geographic Coordinate System and projection for a feature (the feature currently has an Undefined GCS) that I created from a raster file using the Contour tool in Spatial Analyst. (Although the Undefined data lines up with other (defined) features that I have available, as I understand it, it is important to have all features within a map using the same GCS and projection? Is this correct, or can I work across multiple GCS?)

I attempted to set the GCS to GDA_1994, and the projection to Lambert Conformal Conic (these are the GCS and projection of the other data sets that I have available). The data frames GCS is set as above (GDA_1994, Lambert). I also know the approximate location in Decimal Degrees (~ 150.9, -25.0 to 153.9, -28.5 => South East Queensland). When I attempt to use the 'Define Projection' tool (I believe this is the one I should be using, as 'Project' is only used when a GCS is already specified?) the scale and position of the data is completely thrown, with the entire feature fitting within 131.00, 0.00 DD (the measure tool gives the entire data set a size of ~ 3m x 3.5m - pretty much the same values as the true size of the data's range in DD - Coincidence?). Clearly I'm doing something wrong, or I'm missing a step that I have to do before using Define Projection, but I have no idea what..? Do I have to somehow set the datas units, or convert them from the angular units(DD) to linear units(m)?

Any help would be appreciated.


texts:
I'm trying to get map tips to behave as shown below.

Light grey background with very thin white border
A box that adjusts to the content rather than a fixed width of say...300px
A centred title using h1 or h2
Name, company and email in bold with the details all inline with each other (dependant on the longest word in bold).
A small white space
More details in bold with information on the same line not in bold

I'm struggling with aligning [% "Client Name" %] etc to the longest input (shown as "COMPANY NAME") and keeping the bold and normal text on one line.

Current code used:

 ASSET

NAME: [% "First Name" %]   [% "Surname" %]

COMPANY NAME:  [% "Company Na" %]

Progress Report: Transformed into table format

ASSET

NAME:	 [% "First Name" %] [% "Surname" %]

COMPANY NAME:	 [% "Company Na" %]

MANAGEMENT DATE:	[% "Mangement" %]


texts:
Using Qgis 2.0.1, I'm attempting to interpolate a large point dataset (n > 400K) to a raster surface.

I'm using the IDW option in Raster>> Interpolation>> Interpolation, and the resulting output surface has an inexplicable shift 100 metres southwards of the orginal point grid. Can't find the cause of the problem. Has anyone encountered anything like this before?


texts:
I am trying to use GRASS GIS for the first time.
After quite a lot of trouble I managed to create a location and a new mapset, and my shapefile is now displayed in the window "GRASS GIS Map Display" (for some reason only in 3d view).

My next step would be to use the v.to.3d module to add a value to the z-coordinate. I don't seem to get it right, first I tried to use "Search module", and when I found v.to.3d I doubleclicked and filled out the appearing form. That didn't work...
Then I tried writing it manually in the "Command console", but being a total GRASS GIS-beginner I'm not sure how to write? I get, or I think I get, that I need the command 'v.to.3d', an input, output and in some way add a height value.

Can someone please give me a hand? Exactly how can I do this?


texts: Updating shared custom functions / expression code in QGIS

texts:
All the QGIS documentation I've seen on custom functions (used in expressions) says that the code is saved in the user's local expressions directory, e.g. ~/.local/share/QGIS/QGIS3/profiles/default/python/expressions in Linux. The same documentation typically says something like "remember to share your expression files if others need this functionality" (example)

It seems like the only way to share this code is for other users to manually copy it to their own expressions directory, however that seems inadequate. If multiple map developers share a version control system the onus is on them to manually copy any updates to custom functions to their expressions directory. It does not seem possible for expressions to be loaded from an alternate path, e.g. one relative to the .qgs file, and therefore automatically update each time a project is loaded.

I plan to create a Feature Request for loading custom functions from a relative path, but first wanted to check that I have understood the current situation correctly. Do all users of shared custom functions really need to manually copy the code to their expressions directory each time it changes?


texts: Getting higher stream number for lower order stream in ArcGIS Desktop?

texts: The fourth order stream is showing greater stream number then third order stream?


texts: Hillshades created from QGIS Virtual Raster Tables sometimes have unwanted artifacts

texts:
Using QGIS 2.14.3 I tested Virtual Raster Table technology for the first time using two TIF DEMS (NED 1/3 arc-second, projected to EPSG:3338, single band, 32-bit floating point, ~ 5-meter resolution).  I then created a GeoTIFF hillshade from the VRT (using Raster > Terrain Analysis > Hillshade).  I think that the results are spectacular (Brooks Range, northern Alaska):

Happy with the above test, I then used the same procedures to create a VRT for my entire project area (77 NED TIFs), including a hillshade.  Sadly, though, this hillshade (~ 70 GB!) has an unwanted grid pattern:

The same process is applied in both cases, the only difference (that I can think of) is the number of TIFFs in the VRT.  The grid pattern reminds me of improper resampling techniques that I've seen in the past (such as using Nearest Neighbor when Bilinear would have been appropriate), but even still, why did it work so well in the test case?  Any suggestions to avoid this grid pattern?

EDIT: After further experimentation, I have found that this unwanted grid pattern arises when as few as two tiles are used in a VRT.


texts: Layer Syling icons and font size QGIS 3.8.3 64 Bit

texts:
I've installed the 64 bit version of QGIS 3.8.3
I use an external full HD screen. Layer styling fonts and icons are extremely big. I've changed the font and icon size in options, but the problem persists.
Is there a solution to this?

Attached is a picture of the issue.


texts: How to make polygon shapefile containing WGS 84 UTM zones as features?

texts:
I need a shapefile with polygon features. Each feature has to be a unique UTM zone (each zone is 6 degrees) in the northern hemisphere.

I have no idea how to make it - Create fishnet gives a wrong result, and making it by hand is a time consuming operation.

How to automate this, or maybe download shapefile somewhere?


texts:
I have an ASCII grid (Data information in this box link https://app.box.com/s/jpdalfjnyfhqbvrghd8csu3pz5b2fl2j) with the following projection information

Projection    POLAR
Units         METERS
Zunits        NO
Xshift        0.0
Yshift        0.0
Parameters    6378137.0  6356752.314245179
   0  0  0.0 /* longitude of central meridian
  45  0  0.0 /* latitude of true scale
0.0 /* false easting (meters)
0.0 /* false northing (meters)

 PROJCS["Stereographic_North_Pole", GEOGCS["GCS_WGS_1984",
 DATUM["D_unknown",SPHEROID["WGS84",6378137,298.257223563]],
 PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],
 PROJECTION["Stereographic_North_Pole"],
 PARAMETER["standard_parallel_1",45], PARAMETER["central_meridian",0],
 PARAMETER["false_easting",0], PARAMETER["false_northing",0],
 UNIT["Meter",1]]

How can I plot it using the projection information correctly ?


texts: Geoserver URL rewrite is not working in IIS

texts:
I have a website which is running on default port 80 in IIS and it's https, and also in the same server , a geoserver instance is running on port 8080.

I need to redirect all the geoserver requests from my website to geoserver(8080) . I wrote a URL rewrite in IIS , but its not working ,

       <rewrite>
            <rules>
                <rule name="geoserver">
                    <match url="^geoserver/(.*)" />
                    <action type="Rewrite" url="http://localhost:8080/geoserver/{R:1}" />
                </rule>

            </rules>
        </rewrite>

And in geoserver global configuration ,I've put proxy base url as following ,

https://mydomain.com/geoserver/


texts:
Is there a QGIS tool that does the function of http://pro.arcgis.com/en/pro-app/tool-reference/spatial-statistics/grouping-analysis.htm

MY NEED - I have about 15,000 points. I need them clustered in batches of about 700 and polygons created for each batch. There will be a central node that it needs to radiate from.

This was asked back in 2013 in Clustering of multiple raster layers to create sub-regions using QGIS? but there were no replies.

This should look like:

There is almost a PostGIS (and other) solution in the following but this is now ~4y old
Identifying clusters in vector point data using QGIS?

I also had a look at the following but can't find anything that does what I need.
http://opensourcegisblog.blogspot.com.au/2015/03/spatial-analysis-with-qgis-part-i-point.html

I NEED a result like this. Is this essentially the scipy k-means module where the clustering process has several issues?


texts: Move GeoTIFF file to another coordinate of lat and lon

texts:
I have a GeoTIFF file that corresponds to terrain data (SRTM 3-arc data).
I need to do 3 things;

Crop a specific area of the raster and
Literally move the raster to another coordinate of lat and lon.
Save it into new GeoTIFF

For example, if it were a set of islands, cut only one, and move it to another location in the world.

I have QGIS for iMac.


texts: Counting number of times value/string occures in field using ArcGIS field calculator?

texts:
How do I count the number of times a value/string etc. occures in a field (using field calculator)?

In my example I want to count the number of times the letter "H" occures in a field; see attached image.


texts: Predicting the raster from the time series of NDRE using non-linear regression model in R

texts:
Here I have posted the code and got the stacked image of NDRE index by using get_filenames function and then I stacked, there are 4 layers in the stacked image I want to come up with the predicted output. I also has other 2 stacked indices like ndvi and reip. I want to add these indices in the model as well. In the following code I have come up with the linear model like. But I want to perform a non linear model in this for example GAM model.
library(raster)
get_filenames <- function(path,pattern_str="*_aggreagted.gpkg",recursive = F) {
  list_of_filenames <- list.files(path, pattern = glob2rx(pattern_str),
                      full.names = TRUE, recursive = recursive)
  return(list_of_filenames)
}
ndre_16 <-get_filenames("/stacked indices", pattern_str = "stackindicesNDRE.tif")
mmm <- raster::stack(ndre_16)

## create time information
timestamps <- 1:nlayers(mmm)

## slope and intercept
lm_intercept <- calc(mmm, fun = function(x) {
  if (all(is.na(x)))
    return(NA)
  else
    return(coef(lm(x ~ timestamps))[1])
})

lm_slope <- calc(mmm, fun = function(x) {
  if (all(is.na(x)))
    return(NA)
  else
    return(coef(lm(x ~ timestamps))[2])
})

## predict ndvi (y = m*x + t)
ndrepre <-lm_slope * mmm + lm_intercept
plot(ndrepre)


texts: Using Photo Monitoring plugin in a macro in Fiji (ImageJ)

texts:
I am trying to run a macro that will automatically open my 8 band geotiff, stack it has an RGB(false color with Red, Green and NIR) and at the end use the Photo monitoring plugin to create the NDVI. I have this:
run("Bio-Formats Importer","open="+thermal_path +" color_mode=Default rois_import=[ROI manager] specify_range view=[Standard ImageJ] stack_order=Default c_begin=2 c_end=4 c_step=1");

run("Enhance Contrast", "saturated=0.35");

rename("temp_1");

run("Stack to RGB");

rename("for_NDVI");

run("Single image index processing from displayed image", "  select=[NDVI: (NIR-Vis)/(NIR+Vis)] display minimum=-0.50 maximum=0.80 saturation=2 channel=green channel_0=blue select_0=NDVI_VGYRM.lut save");

Everything works up to the use of "Single image index processing from displayed image" tool which gives me multiple Java errors. It works manually as I am simply recording this step using the record macro tool.


texts: How to avoid field duplication when joining two vector layers?

texts:
I have 2 point vector layers wich differ in 2 fields while both of them have the same number of fields (10).
What I want to do is produce a new point vector layer with 12 fields.
When I join the two layers I have as a result a new vector layer which holds all fields from both attribute tables.
Ideally, I would like to join the layers and not keep the duplicates.


texts: Sorting MXDs with specific attributes and exporting to specific directory?

texts:
I'm working on creating a GP tool that simply exports all MXDs in an input folder into an output folder as PDFs. This actually works, but I messed something up with output directories. I coded something incorrectly so where both loops export to different folders, and not to the output folder. Can someone check this and see what I'm doing wrong with my output directories?

This worked well until one of the MXDs has data driven pages enabled, which exports only the first sheet. I believe this is because data driven pages can't be exported using [mapping.ExportToPDF][1], but instead with [datadrivenpages.exportToPDF.][2] To fix this, I test each MXD to see if they have an attribute 'datadrivenpages'. If yes, pass it through the ddp function, else pass it to the normal export function. I messed up my output directories, but I cant figure out how to fix it!

# Broken code here. Working code is below.
# Script takes all .mxd files in a user defined directory and
# exports them to another user defined directory
# This script accounts for map documents with data driven pages
import os, arcpy

# user sets input and output directories
input_dir = arcpy.GetParameterAsText(0)
    #C:\Users\me\Documents\Projects\City Construction\mxd
output_dir = arcpy.GetParameterAsText(1)
    #C:\Users\me\Documents\Projects\City Construction\pdf
arcpy.env.workspace = ws = str(input_dir)

list = arcpy.ListFiles("*.mxd")
mxdnames = []
for item in list:
    print(str(item))
    mxdnames.append(str(item))
print('All Mxds - ' + str(mxdnames))

for mxd in mxdnames:
    input = os.path.join(ws, mxd)
    mxd = arcpy.mapping.MapDocument(input)
    output_name = (os.path.join(output_dir, mxd[:-4]) + ".pdf")
    if hasattr(mxd, "dataDrivenPages"):
        ddp = mxd.dataDrivenPages
        ddp.exportToPDF(output_name, "ALL")
        arcpy.AddMessage('\t' + str(mxd) + str(' exported successfully (Data Driven Pages Enabled).'))
        del mxd
    else:
        arcpy.mapping.ExportToPDF(mxd, output_name)
        arcpy.AddMessage('\t' + str(mxd) + str(' exported successfully.'))
        del mxd

arcpy.AddMessage(('Exports complete!'))
complete = os.startfile(output_dir)
complete

----------

# Solved
# This script prompts a user for an input directory containing MXDs,
# and prompts for an output directory where PDFs are batch exported
# This script exports both normal MXDs and MXDs with data driven pages enabled.

import os, arcpy

# user sets input and output directories
input_dir = arcpy.GetParameterAsText(0)
output_dir = arcpy.GetParameterAsText(1)
arcpy.env.workspace = ws = str(input_dir)
mxdnames = arcpy.ListFiles("*.mxd")

for documentname in mxdnames:
    input = os.path.join(ws, documentname)
    mxd = arcpy.mapping.MapDocument(input)
    output_name = (os.path.join(output_dir, documentname[:-4]) + ".pdf")
    if hasattr(mxd, "dataDrivenPages"):
        ddp = mxd.dataDrivenPages
        ddp.exportToPDF(output_name, "ALL")
        arcpy.AddMessage('\t' + str(mxd) + str(' exported successfully (Data Driven Pages Enabled).'))
        del mxd
    else:
        arcpy.mapping.ExportToPDF(mxd, output_name)
        arcpy.AddMessage('\t' + str(mxd) + str(' exported successfully.'))
        del mxd

arcpy.AddMessage(('Exports complete!'))
complete = os.startfile(output_dir)
complete
````````````````````````````````````
[1]: https://desktop.arcgis.com/en/arcmap/latest/analyze/arcpy-mapping/exporttopdf.htm
[2]: http://desktop.arcgis.com/en/arcmap/10.3/analyze/arcpy-mapping/datadrivenpages-class.htm


texts: Splitting and ordering lines by polygon intersection?

texts:
I am using PostGIS and QGIS 2.2.

I have a shapefile of rivers and a shapefile of state parks and other polygons.  I would like to show the user a simplified version of a stream as a straight line, with the edges of the stream line highlighted where it intersects a polygon.

Here is a picture of my map in QGIS:

Ignore the tributaries - they are well named and easily culled :)

Here is the simplified version I would like:

In order to create the above image, I would imagine I need to split the stream into several lines. In the above example, I would need six.  Blue, Green, Blue, Green, Blue, and Green again.  Each of these segments would have a distance and a type (blue, or green). The total distance would be normalized to unit length, and all the segments would be displayed for the user.

Another solution would be to simply track the Green segments, each with two offsets: a point, and an end point.  In this case, there would only be three segments, and they would show where they started and stopped. This would work well if I wanted multiple types.

It was fun to find the distance of the stream inside of polygons vs outside using ST_Intersect and ST_Length. However, that was simply a ratio of inside vs outside. My solution needs both length and order. I think I need to split the line into several smaller lines - both inside and outside the polygons.

I think my solution needs to satisfy the following properties:

Length of each segment must be preserved
Order of segments must be preserved


texts: Complete Line vector with end points in point shape

texts:
I have 2 different shapes:

lines representing ducts
points representing manholes/poles

Goal:
Fill in 2 different columns start_point and end_point of each line/duct where it touched the points shape.
Each points element has an ID that should be inserted in lines shape, identifying where the duct starts and ends.
Can you point me to the correct direction?


texts: Is there a way to reshape graphics in ArcEngine?

texts:
Our GIS-tool should allow the customer to reshape self drawn (polygon) graphic elements, like one can  reshape features in edit mode with the Editor.

Is it possible to reshape graphic elements?

If not, I think of turning the polygon graphics to features.

Update:

I would accept a simple 'No' :)


texts: Search for all points within a bounding box

texts:
I've already checked all previous similar question. I'm asking again just because I need to be more specific by describing the database structure I am working with.

CREATE TABLE stores (name varchar, geom geometry);
CREATE INDEX geoid ON stores USING gist((geom::geography))

INSERT INTO stores VALUES ('Store name', 'ST_SetSRID(ST_MakePoint(16.123, 51.321), 4326)')
/* where 16.123 is longitude and 51.321 is latitude */

Until now we only perform search-in-radius (circle) queries, with the following sintax

SELECT ST_X(geom) as longitude, ST_Y(geom) as latitude FROM stores WHERE
ST_DWithin("stores"."geom", ST_MakePoint(15, 50)::geography, 10000)
/* where 15 is longitude, 50 is latitude and 10000 is the radius in meters */

As stated in the title, now we also need to perform some query to get all stores within a bounding box, starting from lat_north, lng_east, lat_south, lng_west.

Should I use the ST_DWithin function again? Maybe with the ST_MakeEnvelope. I'm also confused by the use of the geography cast in the "radius" query.


texts: Invalid Geometry, what's wrong with the geometry?

texts:
I am using QGIS (2.10) to make a growing zone map of the U.S. from some data I found online.  I need to reshape to fit into a U.S. map, and I find it's easier if I break the shape into smaller parts, and run an intersect on the smaller part with the map of the U.S. that I have (takes less processing time to do this, then merge the smaller with the larger after).  The problem is that it will not let me break the shapefile, because it is an "Invalid Geometry".  I have brought it into GRASS and have used v.clean, to no avail.

Is there a function that will actually tell you why the shape is an invalid geometry, or can someone tell me how to purge the invalid geometry demon?

I do not see how to attach a file to this post or I would have attached the shapefile.


texts: How to make world file in WGS84 out of OpenLayer's Google Earth layers?

texts:
I would like to make world file out of OpenLayer's Google Earth.
The problem is, that Google Earth layers in Qgis are shown EPSG:3857, but I need just standard WGS84 (ESPG:4326). When I try to set layer units to decimal degrees, Qgis does not change anything. If i force WGS84 projection, than map displayed incorrectly.

Could you please help with this ? I couldn't imagine, that this "map stuff" is so complicated...

Actually, Qgis is not necessary for me... I just need a simple way to get map with world file (WGS84). Are there some small free applications, which could take Google Earth image and make world file for them ?


texts:
I have a dataset which contains location and signal strength data from a field survey for a radio project I'm working on. An example:

lat,      lon,      signal
34.123,  -89.123,  55
34.124,  -89.124,  50
34.125,  -89.125,  36

Of note: the data doesn't fully cover the region it's collected from, because the signal strength is collected at the base station. When the mobile unit can't contact the base, there's no data.

So, the visualization I'm looking for is a sort of combination heatmap and interpolation—the visualization should fade out when there are no points in a given location, like a heatmap, but its color should be determined by the value of the collected points, like an interpolation.

How can I accomplish this?


texts:
I have a question about sensing data time zone from https://scihub.copernicus.eu/dhus/#/home.
In which time zone copernicus provide for us data ?
For example on image we can see Sensing Date Is it CET time zone or what ? Could anyone help me with this question ?


texts: Plot large batch of points without coodinates or geocoding

texts:
I have a list of 50,000 addresses with no Lat/Lon data. I would like to plot these points over a boundaries shapefile. Most geocoding services limit you to 2,500 - 15,000 geocode without paying.

Is there a way to plot these points without geocoding? Perhaps turn the list of addresses into some KML file?


texts: Census.gov is down for now.   I have a friend who is a teaching a GIS class in the San Francisco bay area and only has so many datasets downloaded for specific projects.  Does anyone have a large copy of this area or know of any mirror that the data could be downloaded from?  The class is using QGIS and is going to stall or have to start using data from another source for now to continue.  Any suggestions on other data sets to explore?


texts:
I am very new to using Python...

I have an MXD with all of my layer files in it. They are grouped by category etc. Each layer file references a dataset within my SDE. My SDE has two different connections, Admin and View. Admin allows administration of all feature classes. View is a select only connection.

All of my layers should be connected to feature classes in the View SDE connection. However I have over 70 layers which reference the Admin SDE connection. How do I change this using Python. I do not want to manually change the source SDE connection and then save a new layer.

Would it also be possible to overwrite the layer file using its current file information and location?

I have tried the following bits of code but no luck...

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
for lyr in arcpy.mapping.ListLayers(mxd):
     lyr.replaceDataSource(r"C:\Users\x939062\AppData\Roaming\ESRI\Desktop10.2\ArcCatalog\MMO_View.sde", "SDE_WORKSPACE",,"False")

I have also tried the following, but again no luck...

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
for lyr in arcpy.mapping.ListLayers(mxd):
    lyr.findandreplaceworkspacepath (r"C:\Users\x939062\AppData\Roaming\ESRI\Desktop10.2\ArcCatalog\MMO_Admin.sde", r"C:\Users\x939062\AppData\Roaming\ESRI\Desktop10.2\ArcCatalog\MMO_View.sde")


texts: Could not join CSV to shapefile using PyQGIS 3.16.6

texts:
vlayer = QgsVectorLayer("C:/yyy.shp", "yyy.shp", "ogr")
shp_uri ='C:/yyy.shp'
shp = QgsVectorLayer(shp_uri, 'yyy.shp','ogr')
QgsProject.instance().addMapLayer(shp)
csv_uri="file:///C:/xxx.csv?delimiter=,"
csv = QgsVectorLayer(csv_uri, "Output","delimitedtext")
QgsProject.instance().addMapLayer(csv)
shpField='Piezometer'
csvField='piezometerid'
joinObject = QgsVectorLayerJoinInfo()
joinObject.joinLayerId = csv.id()
joinObject.joinFieldName = csvField
joinObject.targetFieldName = shpField
joinObject.setUsingMemoryCache(True)
joinObject.setJoinLayer(csv)
shp.addJoin(joinObject)

In the joined shapefile, the fields of csv were added, but with NULL values.
I have validated the shpField and csvField.
Also validated the field type for both shpField and csvField.
Is it the problem in my code or the implementation?
If yes, please tell what it is.


texts: Easiest way to combine multiple polygon records into one geometry value?

texts:
I have multiple spatial records in a Postgres table, imported from a shapefile using PostGIS. I am also using QGIS to view & edit the spatial layers.

What is the easiest way to combine the geometry for multiple records into a single geometry value?


texts: How to reduce the size of a NetCDF in R?

texts:
I'm working with a large land cover file (~2.3 GB) saved in NetCDF format (it's the 2018 gridded map from here). My aim is to associate GPS points with this land cover data using R.

The problem is that R cannot allocate a vector of sufficient size to get at the land cover data types once I load it.

Is there a way to bound the data first (I want only the land cover from sub-Saharan Africa) before I try to manipulate it?

Here's my code so far with a failed attempt (based on this tutorial):

#' load any packages
library(ncdf4)

# set path and filename
ncpath <- "C:\\Users\\Me\\Documents\\Science\\GIS\\Global Land Cover\\"
ncname <- "C3S-LC-L4-LCCS-Map-300m-P1Y-2018-v2.1.1"
ncfname <- paste(ncpath, ncname, ".nc", sep="")
dname <- "lccs_class"

# open a netCDF file
ncin <- nc_open(ncfname)
print(ncin)

# get longitude and latitude
lon <- ncvar_get(ncin,"lon")
nlon <- dim(lon)
head(lon)

lat <- ncvar_get(ncin,"lat")
nlat <- dim(lat)
head(lat)

print(c(nlon,nlat))

#' create a bounding box to work with a subset
LonIdx <- which( ncin$dim$lon$vals > 30 & ncin$dim$lon$vals < 32)
LatIdx <- which( ncin$dim$lat$vals > -27 & ncin$dim$lat$vals < -26)
MyVariable <- ncvar_get( ncin, dname)[ LonIdx, LatIdx] # THIS THROWS AN ERROR
# 'Error: cannot allocate vector of size 31.3 Gb'


texts: Label only selected features with PAL PyQgis using an expression

texts:
I'm trying to label only selected features with an expression using PyQgis and the PAL functions, but I don't find any example and the API don't help me a lot (or I don't found it). My Qgis version is 2.14

This is an extract of my code:

from PyQt4.QtGui import QColor
layer = iface.activeLayer()
palyr = QgsPalLayerSettings()
palyr.readFromLayer(layer)
palyr.enabled = True
palyr.bufferDraw= True
palyr.bufferColor= QColor("white")
palyr.bufferSize = 1
palyr.scaleVisibility = True
palyr.scaleMax=2000
palyr.isExpression = True
# code is a fieldname. Only show the last 3 letters
palyr.fieldName = 'right( code ,3 )'
palyr.size = 20
# next 2 lines is for change placement
palyr.placement= QgsPalLayerSettings.OverPoint
palyr.setDataDefinedProperty(QgsPalLayerSettings.Size,True,True,'8','')

# I want to set this expression for label only selected features:  "code not like '%000'"
# MISSING CODE FOR SET DEFINED LABEL USING EXPRESSION

palyr.textColor = QColor("red")
palyr.drawLabels = True
palyr.writeToLayer(layer)
iface.mapCanvas().refresh()

I can run this code without errors, but I don't know how to set expression for defined data features.


texts: how to avoid Intersection of polygons via python code

texts:
i am using QGIS 1.8 and i can set the check box in Setting\Snapping Options "Snapping and Digitizing " dialog to avoid polygon intersections and enable topological editing but i am looking for a way to achieve the same via python code. is there any way to do it with python ?.

Thanks ,
Mael


texts: Calculate traveling time from several points to a same target, for different means of transport

texts:
Provided the following dataframe in which I have several source locations and destination (same destination for all, in this case):
lon_from <- c(2.885882, -2.885882, -2.885521, -2.887446, -2.885882),
lat_from <- c(43.291695, 43.291695, 43.296608, 43.298549, 43.291695),
lon_to <- rep(-2.885882, 5),
lat_to <- rep(43.297457, 5)

df2 <- data.frame(lon_from, lat_from, lon_to, lat_to)

> df2
   lon_from lat_from    lon_to   lat_to
1  2.885882 43.29169 -2.885882 43.29746
2 -2.885882 43.29169 -2.885882 43.29746
3 -2.885521 43.29661 -2.885882 43.29746
4 -2.887446 43.29855 -2.885882 43.29746
5 -2.885882 43.29169 -2.885882 43.29746
>

I would like to calculate how long does it take from each source (c(df$lon, df$lat)) to their target destination (c(df$lon_to, df$lat_to)), according to different means of transport (let's say: car and walking).
I know openroutservice has an R package which has a function, ors_matrix(), which according to official documentation, generates:

One to many, many to many or many to one: ors_matrix() allows you to
obtain aggregated time and distance information between a set of
locations (origins and destinations). Unlike ors_directions() it does
not return detailed route information. But you may still specify the
transportation mode and compute routes which adhere to certain
restrictions, such as avoiding specific road types or object
characteristics.

That seems to be exactly what I am looking for!
If I follow the description, I get pretty interesting results:
library(openrouteservice)
# openrouteservice::ors_api_key("<my_api_key>")

from <- list(c(-2.885882, 43.291695), c(-2.885882, 43.296608))
to <- c(-2.88276, 43.29196)

res <- ors_matrix(from)

> res$durations
       [,1]   [,2]
[1,]   0.00 142.12
[2,] 138.86   0.00

Unfortunately, I haven't been able to tame it for my case scenario. Currently I am facing two problems:

I need to feed it with a list, not a dataframe (although documentation specifies that I could use a dataframe, I haven't been able  to do it).
Distances are calculated from every observation to every other one. I haven't found the way to calculate the many to one option stated in the description.

So, here I am:close to what I am looking for, but stucked. Do you happen to know how to overcome these two problems or any other way to achieve what I am looking for?


texts:
I've been trying to list al file paths in my project using PyQgis (I'm not really acquainted with it); my goal is to translate and clip them all to export them to a folder
I've tried with:

names =    [layer.name() for layer in QgsProject.instance().mapLayers().values()]

But that only returns the names on the layer panel, I want the full path (and to export it as a csv would be great)


texts: SQL Server STRelate equivalent of Oracles SDO_RELATE(...mask=anyinteract)

texts:
In Oracle Spatial, we are running checks on a graph. The Check should provide Nodes, which are marked as endnodes (Attribute "EDGEDEGREE" equals 1) but are in fact located somewhere along an edge. We are using Oracles Relate function to do this with the spatial relation pattern "mask=anyinteract"

In Oracle the Query looks like this:

SELECT 'Node' AS OBJECTTYPE,n.objectid AS OBJECTID,...
  FROM NODE n, EDGE e
  WHERE n.edgedegree=1
  AND n.neighbornode_objectid IS NULL
  AND SDO_RELATE(n.SHAPE,e.SHAPE,'mask=anyinteract')='TRUE'
  AND e.nodefrom_objectid<>n.objectid
  AND e.nodeto_objectid<>n.objectid;

Now we want to do the same on SQL Server (2012 has advanced spatial relation functions now). STRelate uses DE-9IM Sequences to do this. mask=anyinteract corresponds to "non-disjoint" or simply "intersect", but there are 4 different types of intersect. Which DE-9IM Sequence corresponds best to mask=anyintersect for the following T-SQL Query?

SELECT 'Node' AS OBJECTTYPE,n.objectid AS OBJECTID,...
  FROM NODE n, EDGE e
  WHERE n.edgedegree=1
  AND n.neighbornode_objectid IS NULL
  AND n.SHAPE.STRelate(e.SHAPE.MakeValid(),'<?DE-9IM Sequence?>')='TRUE'
  AND e.nodefrom_objectid<>n.objectid
  AND e.nodeto_objectid<>n.objectid;

Could you please explain me these four intersection possibilities, when relating a point on a line?


texts: I tried to clip a USGS GMTED2010 elevation data raster in order to match it with a raster I clipped from the Google Earth satellite. However, when I run the command, the output raster appears in different location with different co-ordinates.