texts: Re-projecting a Raster in R to matching resolution, crs, and origin, but different extent

texts:
I have raster layers in a few different projections that I've loaded into R. I eventually plan on using the "merge" function of the raster package to bring these together as one surface, but to do that, all of the rasters need to have the same crs, resolution, and origin. Using projectRaster with "alignOnly=TRUE, says it will use "to or other parameters only to align the output (i.e. same origin and resolution), but use the projected extent from from," which sounds like what I need. However, this fails and gives me the warning "Warning in projectRaster(x, crs = projection(y)) : 'from' has no cell values." Does anyone have any suggestions on how to accomplish this? Is this a bug, or am I interpreting the documentation incorrectly?

Below is a small reproducible example of trying to reproject one raster with the GRS80 ellipsoid and 1.5m resolution to align with another raster of a different extent that uses the WGS84 ellipsoid and has 4m resolution.

#Create a raster with GRS80 crs and 1.5m resolution
r1<- raster(xmn=731458, xmx = 731458+100*1.5,  ymn=3116234, ymx= 3116234+100*1.5,
nrow=100, ncol=100, vals=sample(1:50, size = 100*100, replace = TRUE), res=1.5,
crs= "+proj=utm +zone=16 +ellps=GRS80 +towgs84=0,0,0,-0,-0,-0,0 +units=m +no_defs ")

#Create a raster with WGS84 crs and 4m resolution
r2<- raster(xmn=740000, xmx = 740000+100*4,  ymn=3000000, ymx= 3000000+100*4,
nrow=100, ncol=100, vals=sample(1:50, size = 100*100, replace = TRUE),
res=4, crs= "+proj=utm +zone=16 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")

#Reproject r1 to WGS with 4m resolution
r1_res4_WGS84<- projectRaster(from = r1, crs = crs(r2), res= res(r2))

#Origins do not match after reprojecting
origin(r1_res4_WGS84) #1.999999e+00 9.094551e-05
origin(r2) #0 0

#Try reprojecting with alignOnly=TRUE
r1_res4_WGS84_align<- projectRaster(from = r1, to = r2, alignOnly = TRUE)
        # Warning in projectRaster(x, crs = projection(y)) :
        #   'from' has no cell values


texts: What are the pros and cons of running on a windows server vs. a linux server?

texts:
We're in the process of trying to decide where to host our GeoServer installation in production.

Are there any major pros or cons of hosting on windows vs. linux?


texts:
I am having a strange problem--currently using QGIS 3.03. I cannot zoom to any of my point layers (can't use zoom to selected, zoom to layer, etc.)  But I CAN zoom to my polygon layers.  Both are shapefiles, both are WGS84 unprojected (EPSG4325). Does anyone know what might be going on?

I saw a related question that said to disable a plugin called Tile Map Scale, but I have no such plugin installed.  Just as an experiment, I disabled and uninstalled all plugins, and am still having the same issue.

Are there perhaps any other zoom or scale settings I might be overlooking?


texts: Delete records in feature class that have same ID in CSV table using arcpy cursors

texts:
So I have a csv table with records that have an ID.  I want to find the same IDs in a feature class/shapefile and delete those records in the feature class/shapefile using Python.

I've messed with add.join and relates without a good result so I figured I'd move to cursors.  I've setup two search cursors to test the result but I don't think the logic is quite right so I'm posing my question here.

Here's a code snippet.  The second cursor will eventually change to Update.Cursor using delete row but not until I get the result correct first.

import arcpy

tc = "d:\\temp\\tc.shp"
changes = "d:\\temp\\changes.csv"

with arcpy.da.SearchCursor(changes,"RECORDNUM_AFTER") as changecur:
    for row in changecur:
        changeIDs=row

tcur = arcpy.SearchCursor(tc, fields="RECORDNUM")
for row in tcur:
    if row in changeIDs:
        print row.getValue("RECORDNUM")


texts: How do I develop my GIS programming skills?

texts: I would like to develop my GIS programming skills, where do I start? People say learn C++ or Python but where can I learn this in a geographic context? Tutorials or anything would be very helpful, as would any information on languages/programmes to use.


texts:
I am attempting to create a single set of all the attributes read from a field on a feature class.  This is my code so far:

import arcpy
workspace = arcpy.mapping.MapDocument("CURRENT")
cursor = arcpy.da.SearchCursor("FireStations",['CITY'])
s1 = set()

for row in cursor:
    s1.add(row)
    print(s1)

The issue I am having is that it seems to create a new set with all the attributes within the field for every new attribute.  So instead of one single set of all the attributes I get one set with all the attributes for every row in the specified field. I know this is most likely an issue with my for loop, anyone have any thoughts on fixing this issue?


texts: Republishing feature layers in ArcGIS Online but getting ERROR 001270?

texts:
Trying to republish feature layers in ArcGIS online but getting an ERROR 001270. Basically, we were having issues exporting/downloading a feature layer to FGDB from Arc GIS Online. ESRI support stepped in because nothing seemed to be working and the data had information and attachments that we couldn't just get rid of. They sent us the FGDB and I copied it and edited it with the hopes of republishing it. I checked/repaired the geometry and tried basically every work around I could find online regarding error 001270. The feature is a point file that has a relationship class with the attachments table but now you just see the table and feature class as separate objects in AGO, not the usual feature with files attached, like before.

I might just have to go back to ESRI support, but I already spent 2 hours of their time getting the FGDB exported and downloaded.


texts: Adding feature class into ArcGIS 10.3 mxd using python script

texts:
I am trying this simple task to add a feature class to an mxd (ArcGIS 10.3) but getting few errors. Here are two options I tried.

1) Preferred option - adding feature class itself into a active dataframe of the current
mxd. I used the script below:

import arcpy

mxd=arcpy.mapping.MapDocument("CURRENT")

DF=mxd.activeDataFrame

arcpy.mapping.AddLayer(DF,"C:\Work\GP1500\Tests\Input.gdb\Poly1","TOP")

The addLayer line generates the error below
......
assert isinstance(add_layer, Layer)
AssertionError

2) Create a layer from the Feature class and add the layer using the script *****

import arcpy

mxd=arcpy.mapping.MapDocument("CURRENT")

DF=mxd.activeDataFrame

FC="C:\Test.gdb\Test_Poly1"

arcpy.MakeFeatureLayer_management(FC,"TestPolyLayer")

inputlayer =arcpy.mapping.layer(TestPolyLayer)

arcpy.mapping.AddLayer(DF,inputlayer,"TOP")

******It generates the error while executing the second last line i.e. to define input layer with following message
AttributeError: 'module' object has no attribute 'layer'

Note: the second option i used is after I saw Nick's reply in How to add feature class to MXD with ArcPy (Python)?


texts: Leaflet tiles are not loading for multiple Bootstrap tabs

texts:
I am making a web app for different areas using Leaflet library. I have different Bootstrap tabs, each tab opens a new map of a new area.

The problem I am facing is that the tiles on the first tab load but the tiles on the second tab do not. Where am I wrong?

This is how I am loading my maps for now...

/* MAP LAYERS */
var map = L.map('mapid',{
        zoomControl: true,
}).setView([24.2, 67.8], 9);

/* LOADING AND ADDING BASE MAPS */
var baseMap = L.tileLayer('http://{s}.basemaps.cartocdn.com/light_all/{z}/{x}/{y}.png', {
                    attribution: '&copy; CARTO contributors',
}).addTo(map);

/* MAP LAYERS */
var map1 = L.map('map_id',{
        zoomControl: true,
}).setView([22.2, 63.8], 9);

/* LOADING AND ADDING BASE MAPS */
var baseMap = L.tileLayer('http://{s}.basemaps.cartocdn.com/light_all/{z}/{x}/{y}.png', {
                    attribution: '&copy; CARTO contributors',
}).addTo(map1);

My tabs:

                        I

                        J

The tab linkage:


texts: How to match UTM squares between Excel and Shapefile?

texts:
I have an excel file with locations displayed in UTM 10kmx10km squares. I also have a shapefile with those same UTM squares and their coordinates, however I cannot seem to figure how to project the excel UTM squares.

I have tried to join both of the attribute tables but with no success. Am I doing something wrong?


texts:
I have the following coordinates in meters for UTM 34S / Datum WGS-84:

Lat.: 311965;
Long.: 8695153;

How do I convert these coordinates into decimal degrees or DMS?


texts: Use PHP form to insert geometry from WKT into SQL Server 2017

texts:
I have a PHP 7.3 form that asks the user for a Name, a WKT and an SRID.  I would like to upload the WKT and SRID to SQL Server 2017 so that it creates a geometric object in Upload_WKT_Test using STGeomFromText:

<?php
    if ($_SERVER['REQUEST_METHOD']=="POST") {
        $wkt = $_POST['wkt'];
        $srid = $_POST['srid'];
        $name = $_POST['name'];
        try {
            $wktQuoted = $pdo->quote ($wkt);
            //$wktQuoted = "'$wkt'";
            $sql = "INSERT INTO Upload_WKT_Test (Name, GeomCol1) VALUES (:name, :wktGeom)";
            $wktGeom1 = "geometry::STGeomFromText(";
            $wktGeom = $wktGeom1."".$wktQuoted.", ".$srid.")";
            //echo $wktGeom."

".$name."

";
            $stmnt = $pdo->prepare($sql);
            $theData = [':name'=>$name, ':wktGeom'=>$wktGeom];
            $stmnt->execute($theData);
        } catch(PDOException $e) {
            echo "Error: ".$e->getMessage();
        }
    } else {
        $wkt="";
        $alignment="";
        $srid="";
    }
?>

My init.php:

<?php
    ob_start();
    session_start();
     try {
        $pdo = new PDO( "sqlsrv:Server=localhost\SQLEXPRESS;Database=devdb", "", "");
        $pdo->setAttribute( PDO::ATTR_ERRMODE, PDO::ERRMODE_EXCEPTION );
        $pdo->setAttribute( PDO::ATTR_DEFAULT_FETCH_MODE, PDO::FETCH_ASSOC );
        //$pdo->setAttribute( PDO::ATTR_EMULATE_PREPARES, true );
    }

    catch( PDOException $e ) {
        //die( "Error connecting to SQL Server" );
        //die(print_r($stmnt->errorInfo(), true));
        echo "Error: ".$e->getMessage();
    }
    $root_directory = "testwkt";
    $from_email = "admin@somewhere.com";
    $reply_email = "admin@somewhere.com";
    include "php_functions.php";
?>

I understand that the preferred method to accept user input via forms is to use a parameterized query to prevent SQL injection.  I believe the source of the error message has to do with the quotes around the WKT are not making it into the INSERT INTO statement.  Is it possible to surround text values from an input form with quotes (WKT) and use it to build a geometry object?

Example WKT entered into form POINT(100 10)

Example SRID 0

Example Name Test

Geometry from WKT (quotes are required):

geometry::STPointFromText('POINT (100 10)', 0)

Resulting error message:

SQLSTATE[42000]: [Microsoft][ODBC Driver 17 for SQL Server][SQL Server]
A .NET Framework error occurred during execution of user-defined routine or
aggregate "geometry": System.FormatException: 24114: The label
geometry::STGeomFrom in the input well-known text (WKT) is not valid.
Valid labels are POINT, LINESTRING, POLYGON, MULTIPOINT, MULTILINESTRING,
MULTIPOLYGON, GEOMETRYCOLLECTION, CIRCULARSTRING, COMPOUNDCURVE,
CURVEPOLYGON and FULLGLOBE (geography Data Type only).
System.FormatException: at
Microsoft.SqlServer.Types.OpenGisTypes.ParseLabel(String input) at
Microsoft.SqlServer.Types.WellKnownTextReader.ParseTaggedText(OpenGisType type) at
Microsoft.SqlServer.Types.WellKnownTextReader.Read(OpenGisType type, Int32 srid) at
Microsoft.SqlServer.Types.SqlGeometry.GeometryFromText(OpenGisType type, SqlChars text, Int32 srid) at
Microsoft.SqlServer.Types.SqlGeometry.Parse(SqlString s) .

Upload_WKT_Test table:

CREATE TABLE dbo.Upload_WKT_Test
    ( id int IDENTITY (1,1),
    Name varchar(50),
    GeomCol1 geometry );
GO

WKT - Well Known Text - a way to represent a geometric object (point, line, polygon for example) in text format.  More info here.

SRID - Spatial Reference System Identifier - an integer that represents a coordinate system.  More info here.

More info about working with Geometry Instances in SQL Server and Azure SQL Database can be found here.

More info about STGeomFromText


texts:
I made a TIN interpolation of some height points but the resulting raster appears with some weird lines (like in the pic)

What can i do to avoid this result?

Here are the input parameters. Sorry that my QGIS is set to Italian.


texts: Changing name layer from geoprocessing output

texts:
How can I change the name for geoprocessing output map layers, for example:
out_diss = 'ogr:dbname=\'C:/GIS/data/testdata/output.gpkg\' table=\"ak_riks_dissolved\" (geom)'

dissolved = processing.run("native:dissolve",
    {'INPUT':'C:/GIS/data/testdata/ak_riks.shp','FIELD':[],
    'OUTPUT':out_diss})['OUTPUT']

dissolved.setTitle('dissolved') #No difference
QgsProject.instance().addMapLayer(dissolved)

Instead of the default very long name:


texts: Find nearest point given coordinates from another table

texts:
So I have two tables, one table is called stores which has the name of the store as well as the geom column that contains the coordinates for the stores. The other table is customers which has the coordinates of the customer's address.

I'm trying to create a new column in the customer table with the nearest store given the customer's coordinates. I've tried something like this but it doesn't work:

UPDATE customers SET storename =
    SELECT storename
    FROM stores
    ORDER BY geom <-> ST_GeometryFromText('POINT(customers.geom)',4326) LIMIT 1
);


texts: How to generate Digital Surface Model (DSM) from aerial image?

texts:
Can I use Erdas LPS 2011 for automated DSM generation with low resolution (5m to 15m) by one aerial image?

I don't have a stereo pair images for the city.


texts: Unable to restore ESRI reference in VS2017

texts: I have project which was developed on ASP.Net 4.5 and ArcGIS Runtime SDK for .Net 10.4. Now I have installed ArcGIS Runtime 10.4 and change framework from 4.5 to 4.6.1 and also tried to install nuget packages but still not able to restore references .I followed instruction mentioned in ArcGIS SDK. How can I restore references?


texts: Converting JPEG2000/MrSid files to GeoTIFF using GDAL

texts:
This question is related to one previously asked: Batch convert jpeg2000 to GeoTIFF?

We certainly need to batch convert some JPEG2000 files to GeoTIFF; the problem is, it doesn't seem to work. We don't have ready access to ArcGIS, so the GDAL solution discussed previously is what we need to work with -- but here's the rub:

The JPEG2000 files in question originated in the MrSID format, and apparently have supporting metadata that make them work properly. When we convert a file with gdal_translate, there are two problems.

First, a GeoTIFF converted in this way doesn't seem to be properly optimized for anything. Loading the file as a layer in QGIS (to test if the output is a valid raster), it takes lengthy minutes to render -- versus a GeoTIFF exported by the LizardTech GeoViewer, specifically designed for this sort of thing, which takes only a second or two to render (even though both files are about 4 GB in size). It would be ideal to be able to batch-convert the files using GeoViewer but there does not seem to be a way to do that.

Second, the GeoTIFF converted by gdal_translate is not able to be served at all by GeoServer. It allows the creation of a layer using the image but then serves only a blank layer. It is unclear why this is happening. Presumably whatever MrSID metadata is present and necessary in the original JPEG2000 is not being carried over to the GeoTIFF.

Any advice on why this conversion is not working, or how to do it correctly?


texts: How to highlight a selected feature in qgis plugin using python

texts:
I am working on a plugin that would highlight selected features by changing its color. Can anyone help me regarding this? The plugin searches for required features and once found, I want that those features are highlighted with color different from the layer color. How to do this?

Here is the snippet:

def processResults(self, results):
    if self.layer is None:
        return
    if self.selectBox.isChecked():
        self.layer.setSelectedFeatures(results)
        if len(results) == 0:
            return

        if self.panBox.isEnabled() and self.panBox.isChecked():
            canvas = self.iface.mapCanvas()
            rect = canvas.mapRenderer().layerExtentToOutputExtent(self.layer, self.layer.boundingBoxOfSelected())
            if rect is not None:
                if self.scaleBox.isEnabled() and self.scaleBox.isChecked():
                    rect.scale(1.5)
                    canvas.setExtent(rect)
                else:
                    canvas.setExtent(QgsRectangle(rect.center(), rect.center()))
            canvas.refresh()

This code just scales the canvas 1.5 times, what I want is not zooming in it, but just highlight it on the same scale by either changing color or any other means. How to do it?


texts: Minimizing bandwidth while refreshing an ArcGIS dynamic map layer

texts:
I've been tasked with creating a web app where multiple users can draw on a map and have that reflected to other users in real time.  Persistence of the drawing if users leave and come back is important.  The server will be ArcGIS Server and the client will be done with the JavaScript API.

I'm pretty sure the best way to do this is with a dynamic map service layer set with a short refresh time, but I'm worried about bandwidth usage.  Are there any strategies for minimizing the volume of data sent such as having the server only send what's changed after the first query?


texts: Layer to KML data interoperability error

texts:
I have a 2D raster dataset that I am trying to convert to KML (using the Layer to KML tool in ArcGIS desktop) for display in Esri Earth, but no matter what I try as the inputs, I keep getting a data interoperability error.  I've tried using a raster layer both within a file gdb and outside of a file gdb thinking this could be causing the issue, nope.  I'm using the defaults for the optional parameters:

Any idea why this isn't working??

code

<?xml version="1.0" encoding="UTF-8"?> <metadata xml:lang="en"><Esri><CreaDate>20170612</CreaDate><CreaTime>17164000</CreaTime><ArcGISFormat>1.0</ArcGISFormat><ArcGISstyle>ISO 19139 Metadata Implementation Specification</ArcGISstyle><SyncOnce>TRUE</SyncOnce><DataProperties><lineage><Process ToolSource="c:\program files (x86)\arcgis\desktop10.4\ArcToolbox\Toolboxes\Conversion Tools.tbx\LayerToKML" Date="20170612" Time="171640">LayerToKML Inundation_ft E:\MP\Inundation_ft_LayerToKML.kmz 0 false # 1024 96 CLAMPED_TO_GROUND</Process></lineage></DataProperties></Esri></metadata>


texts: In Postgis 1.5, how can I split lines in points while retaining the id?

texts:
Hypothesis:
Having 2 linestring tables:
pipelines (id integer, street varchar(50), shape geometry)
streets (id integer, name varchar(50), shape geometry)

Each pipeline passes/intersects several streets.
Need
I need to:

Fill the street field in the pipelines automatically and
I need to know to which street does the pipeline "belong the most".

I thought about densifiying the vertices of the each pipeline at a maximum 5 meter distance from each other, and do the analysis on them and then statistically find out the answer.
How In Postgis can I create a table that presents the geometry for each point of the pipelines along with their corresponding id?


texts: Setting PostgreSQL privileges for ArcGIS?

texts:
I'm migrating from SQL Server and setting up my first PostgreSQL server for use with ArcGIS.  I've followed esri tutorials and related how-tos.  I'll be using a versioned editing environment, as I will want numerous users editing simultaneously.  I'm attempting to set up users/roles following the Esri documentation:

I've created an "editor" role for editing, using pgAdmin.  I've also created an "owner" user for creating versions, using the Create Database User geoprocessing tool.  The owner has an "owner" schema in the database, and I've granted USAGE on the schema to the editor role.  Sufficient privileges are granted to the "editor" role to allow user in the role to edit data.

At this point, the owner can create datasets, feature classes, etc, while the editor can only edit current objects within the dataset.  I've noticed the editor can also create new feature datasets within the geodatabase, but is unable to create new objects in the new feature dataset.

Is there a privilege I can revoke from the editor to disable it from being able to create feature datasets?


texts:
I would like to produce an atlas of 5 blocks for 16 variables.
So for each atlas I display a layer who represent a variable.
For each map in atlas I'm able de set a title with display block name with something like that :

'NDVI ' + "block_name_field"

But I have 16 variables and I don't want to handwrite it. Is code like that possible ?

Active_layer_name + "block_name_field"


texts:
The polygon fill for 'Urban parcels' in the following .gif animation looks really cool.

http://changelog.qgis.org/media/images/entries/5817785cf3d18db78669c2eb2c2be0cf606c3783.gif

(See also: http://changelog.qgis.org/en/qgis/version/2.16.0/)

Does anybody know how to create a polygon fill like this?


texts: How to report issues with QGIS documentation?

texts:
I stumbled across a gotach while installing mod_fastcgi on MacOs 10.8 (mountain Lion) -- there is a bug in 10.8  in the apache build env which breaks the install instructions here: http://hub.qgis.org/projects/quantum-gis/wiki/QGIS_Server_Tutorial

The fix is to put a symlink

 sudo ln -s /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/ /Applications/Xcode.app/Contents/Developer/Toolchains/OSX10.8.xctoolchain

as documented here:

https://apple.stackexchange.com/questions/58186/how-to-compile-mod-wsgi-mod-fastcgi-etc-on-mountain-lion-by-fixing-apxserror

It would be good to get this noted in the docs...

Should I open standard bug report or is there a better way to report 'problems' with a document?


texts:
I am looking for a function with despair in Google Earth Engine that would give me the direction from point 1 to point 2. The goal is to get the same value that Google Earth Pro gives during a measurement (heading).

Does this function exist, how could I get this measurement?


texts: How to snap polygons within a large polygon hole

texts:
I'm trying to perform a snap operation to fill the white holes and remove the intersections made by a previous shapely simplify operation. Here's one example to the occurrence:

There is the code I am running:

import json
from pathlib import Path

from shapely import geometry

POLYS_ROOT_LST = json.load(open("geometries_list.json", "r"))

def create_feature_collection():
    geojson_features = []
    last_shapes = []
    next_shapes = []
    for idx, shapes in enumerate(POLYS_ROOT_LST):
        for shp in shapes:
            shp = geometry.shape(shp)
            medium_point_dist = (shp.length / len(shp.exterior.coords))
            simplify_tolerance = medium_point_dist * 1.2
            shp = shp.simplify(simplify_tolerance)

            # for last_shp in last_shapes:
            #     if shp.within(last_shp):
            #         shp = snap(shp, last_shp, 0.01)
            #         break

            geojson_features.append({
                "type": "Feature",
                "geometry": geometry.mapping(shp),
            })

            next_shapes.append(shp)

        last_shapes = next_shapes
        next_shapes = []

    return {
        "type": "FeatureCollection",
        "name": "Teste",
        "features": geojson_features
    }

if __name__ == '__main__':
    output = Path(__file__).parent / f"test.geojson"
    output.write_text(json.dumps(create_feature_collection()))

To download the geometries here is a link:
https://drive.google.com/file/d/1GOp2v2SqRVG1t8D1G38GC-nO5TJlBKM3/view?usp=sharing

I am using QGIS for the GeoJSON view.


texts: How do I search for features within another feature using WFS OGC

texts:
I understand I can write a GML query like this...

<wfs:GetFeature xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
  xsi:schemaLocation="http://www.opengis.net/wfs http://giswebservices.massgis.state.ma.us/geoserver/schemas/wfs/1.0.0/WFS-basic.xsd"
  xmlns:gml="http://www.opengis.net/gml" xmlns:wfs="http://www.opengis.net/wfs"
  xmlns:ogc="http://www.opengis.net/ogc" service="WFS" version="1.0.0">
  <wfs:Query typeName="massgis:GISDATA.SCHOOLS_PT"
    xmlns:massgis="http://massgis.state.ma.us/featuretype">
    <ogc:PropertyName>NAME</ogc:PropertyName>
    <ogc:PropertyName>ADDRESS</ogc:PropertyName>
    <ogc:Filter>
      <ogc:Within>
        <ogc:PropertyName>SHAPE</ogc:PropertyName>
        <gml:Polygon srsName="http://www.opengis.net/gml/srs/epsg.xml#4326">
          <gml:outerBoundaryIs>
            <gml:LinearRing>
              <gml:coordinates>-71.591935,42.063582 -71.564217,42.063968 -71.563184,42.014226 -71.591207,42.013596 -71.591935,42.063582</gml:coordinates>
            </gml:LinearRing>
          </gml:outerBoundaryIs>
        </gml:Polygon>
      </ogc:Within>
    </ogc:Filter>
  </wfs:Query>
</wfs:GetFeature>

This is fine when you know the geometry ahead of time.  If I have a FeatureId of a feature on another layer, is it possible to include this FeatureId somehow and let the GeoServer work out what the geometry is?  I've spent ages trying to find examples and found nothing like this.  I wonder if perhaps it just isn't possible and I will have to query the server for the geometry, parse it and then send it back to the server.  This seems very longwinded.

Thoughts?


texts: How to create point labels from midpoints of road polylines in QGIS?

texts:
I have a road polyline layer (singlepart ie has several features per road, grouped by an id field), and i want to (after merging the parts) find the mid point of the resulting lines. This is for roads geolocation.

The problem is that the QGIS' dissolve function only works on polygons.

I've tried:
a) mmqgis plugin "create label layer" (handles multiparts nicely)
b) buffer road lines to 1m then dissolve by roadid, then ftools centroids

Both of which give a 'sort of' result, but the main problem is that the point is not actually on the road line.  Is there a way to either force polygon centroids to be within the buffer polygon, or alternatively, some way to actually join the polylines and find mid points directly.

Note that while there are (a small number of) 'Y' shaped lines, most of the multipart lines start and end at the same place, (maybe gridifyed to 1cm).


texts: Raster data does not clip to the vector polygon boundary

texts: I am trying to clip a raster to vector polygon. I am using "Extract by Mask" tool and my output raster layer is, what seems like, randomly cut off in the middle. My processing extent is set appropriately, my coordinate system is the same for both layers, and I have about 10 other raster layers that I did the same thing for and they worked perfectly. I have no idea how to diagnose this problem. Any ideas?


texts: Polygon layer in shapefile not showing which may be projection issue

texts:
It's been a long time since I used GIS software!
I would like to have the symbolised polygons of all Soil Orders on the global map from FAO based on the attribute column "SNAME"
http://www.fao.org/fileadmin/templates/nr/images/resources/images/wsrll.zip
I add the shapefile wsrII to QGIS and I can see the lines and dots, but the polygon features (presumably wsrII_PAL) are not showing up! I assume its a projection issue but I can't work it out.
How do I do this?


texts: Width of a dual-line based on values out of the attribute table in ArcGIS

texts:
Is it possible to symbolize the width of a dual-line automatically by a value out of the attribute table in ArcGIS for Desktop 10.4.1?

In my case, i have hundreds of different entries with the width of a sewer in the attribute table field and wanted to have this scaled distance between the two lines of the dual-line.


texts: Merge GeoNetwork metadata with GeoServer layers

texts: We have a GeoNetwork metadata catalog already running with hundreds of entries and would like to merge it with GeoServer so that the layers in GeoServer have the metadata from GeoNetwork. Is that possible?


texts: Specifying values in stretched raster histograms of ArcMap?

texts:
I am trying to use the stretched color option with a topography raster. I have a color bar that I painstakingly make in arc and I am basically trying to set values of the topography equal to certain colors using the histogram. The problem (I think) is the histogram doesn't have a break where I need it - it has a -40 then jumps to 19. I want a break at 0 so sea level will align with my shift from blue to green in my color bar.

Is there a way to set these values?
Or is there just a better way to do this altogether?

I would use classified color, however it doesn't seem to have a way of ramping the color between the bounds. Its just a single color for a range, and with only a maximum of 32 intervals it doesn't look very good.

Here is a picture of my best result using stretched. You can see the areas on the coast that are green.


texts: Selecting objects based on multiple attributes QGIS

texts:
I would like to select objects from the attribute table based on multiple attributes.

I want to select all ojects that have 'ja' in the column "BLAUWBORST", all objects that have 'ja' in the column "BRUINE_KIE", objects that have 'ja' in both of them and objects that have 'ja' in the column 'veertien' or 'vijftien' or...


texts: How do I set labels to not be one over the other?

texts:
In certain zoom I see many city labels "written" on each other. I know there is a feature that set to auto hide labels that override each other ( and I'm not talking about rendering from specific scale...)
Where can I find it in QGIS?


texts: How to change properties of an item in a legend in ArcMap 10.2 using python

texts:
For many of my projects I use a single mxd to create multiple maps using a python script to switch on/off various raster layers before exporting to pdf. Before beginning the exporting process I currently have to open the legend properties window, go to the items tab, individually double click on each raster layer (often 50+ raster grids) and untick the "Show Layer Name" checkbox for each in order to display it correctly in the legend. As this is a highly manual process I almost always accidentally miss doing this for at least one raster in the list and have to spend time working out which layers to switch on/off before manually exporting.

Is there a way to automate removing the layer name line of each raster in the legend using python? i.e. leaving just the "Heading" and "Labels" boxes ticked.

Screenshot showing the location of the layer name I'm looking to remove and the tickbox to do so below. See red/yellow arrows.

http://gyazo.com/abe72eb654f08845e6f2cfe087ebc613


texts:
I am working on a ArcGIS Python Add-ins which need user login to get data. Rather than exposing user's password in a combobox forever, we need to have a password input box which only show * while user typing.

I tried to replace the "value" of the combobox with stars in script but it only replace all of them after enter. What we need is to replace the character whenever the user type.


texts: GEE: How to transform masked values (NODATA) to -9999?

texts:
How can I transform NoData (masked) values to a value? When the image is mapped, there are a lot of values that are masked and I would like to get a value (-9999) on that masked pixels.

var region = geometry;

// Collect data and filter by total dates
var coll = ee.ImageCollection ("MODIS/006/MOD16A2")
.filterDate('2000-01-01', '2014-12-31')
.filterBounds(region)
.filter(ee.Filter.calendarRange(1,1,'month'));

print (coll)
var modisET = ee.Image('MODIS/006/MOD16A2/2001_01_01')
.select("ET");

var multiply = modisET.multiply(0.1)

var divide = multiply.divide(8)

// Normalize the image and add it to the map.
var rescaled = divide.unitScale(0,1);
var visParams = {min: -1, max: 1};

var reprojected = divide
    .unitScale(0, 1)
    .reproject('EPSG:4326', null, 500);

    Map.addLayer(reprojected)

// Export the image, specifying scale and region.
Export.image.toDrive({
  image: reprojected,
  description: '2001_01_01',
  scale: 500,
  region: table
});


texts: Using "Merge" with in_memory workspace with Arcpy in ArcMap 10.4

texts:
I am attempting to merge three tables, all within the in_memory workspace. The merge appears to work fine, but when I try to call the merged table, it cannot be found. I have isolated the issues to this very basic code:

pt = 'eo_sourcept'
ln = 'eo_sourceln'
py = 'eo_sourcepy'

point = arcpy.TableToTable_conversion(pt, "in_memory", "point")
line = arcpy.TableToTable_conversion(ln, "in_memory", "line")
poly = arcpy.TableToTable_conversion(py, "in_memory", "poly")

sourceFeatures = arcpy.Merge_management([point, line, poly], os.path.join("in_memory", "sourceFeatures"))

arcpy.MakeTableView_management(sourceFeatures, "sourceView")

I get the error:

ExecuteError: ERROR 000732: Input Table: Dataset in_memory\sourceFeatures does not exist or is not supported

I am getting this error for any tool I try to run while using the in_memory merged dataset. If I set the path to a local/permanent drive, the merged table can successfully be called.

The merge tool is the only tool that I have come across that I have had problems with the in_memory output not being able to be found. Any ideas?


texts: Derivation of initial bearing calculation

texts:
NOTE: I was recommended by another user that this would perhaps fit better in this stack, initial post link.
I needed to make calculations of the initial bearing between two coordinates described in latitude/longitude. I have found one formulation of this calculation as

Note, I tried formatting the equation with LaTeX but I couldn't see how to do that on stack, using $$ didn't seem to do anything.
This calculation can be found on various sites like link 1 2 3 but I have not been able to find a derivation or a source for the calculation, as such I am looking for a derivation or an explanation of the derivation.


texts: Combine/Converting many parallel lines to one single line using ModelBuilder of ArcMap

texts:
I am working with railway generalization by using ArcGIS ModelBuilder and my solution needs to involve that software. I need to "Combine/Convert" many parallel lines together into one wherever they are close to each other.
Beside use Merge Divided Road tool, does anybody know how to "Merge" many lines to one?
The left data should become the right.


texts: Measure Snapping with ArcGIS Javascript API

texts: I would like to know if the ArcGIS Javascript API has built in snapping functionality for the measure tool.  I am trying to create a measure tool which only measures distances along the edge of roads.  Is there anything built-in or would I need to create this functionality by linking into ArcObjects?


texts: Is it a good idea to use GRASS Algorithms for Geoprocessing in a Server Environment?

texts:
I´m designing a web-service to run in a server environment that will do some geoprocessing actions to datasets uploaded by users and send the results back.

I´m programming in Python, and so far I´ve been using different GRASS algorithms to do the geoprocessing. Which means, I´m using a lot of code that looks like this:

grass.run_command("v.buffer", input="a", output="b", type="line", distance=c)

Does anyone have experience in using GRASS-algoriths in a server environment? Is it a good idea to use GRASS or should I better be using other libraries (e.g. because they are faster or more stable?). It has to be OpenSource, since we don´t have the budget for ArcGIS Server, etc.


texts: How to export/copy a shapefile from a hosted feature layer from ArcGIS online using ArcPy (Python 3.x) and ArcGIS Pro

texts:
I am rewriting a 2.x script for compatibility with ArcGIS Pro and Python 3.x. My question is, how do I export a shapefile from a hosted feature layer from ArcGIS online using ArcPy (Python 3.x) and ArcGIS Pro?

The working python 2.x script written for ArcGIS Desktop does it this way:

mxd = arcpy.mapping.MapDocument("path_to_mxd_containing_online_layer.mxd")
df = arcpy.mapping.ListDataFrames(mxd)[0]

### exporting the data to a local TEMP layer
inspector_layer = arcpy.mapping.ListLayers(mxd, "name_of_layer_in_mxd", df)[0]
for lyr in arcpy.mapping.ListLayers(mxd, '', df):
    if lyr.isFeatureLayer == True:
        d = lyr.name.encode('ascii', 'ignore')
        file_name = d + ".shp"
        dest1 = os.path.join(wkspce, file_name)
        try:
            arcpy.CopyFeatures_management(lyr, dest1)
        except:
            pass

I was thinking of creating an aprx project file that contained a map with the layer in it, as that seems the closest to the original script. This is one way of writing it:

aprx = arcpy.mp.ArcGISProject("G:\\Code Enforcement\\ACD Admin Operations\\GIS\\Operations\\Mapping\\Automation\\Auto_District_Mapping_PRO\\Auto_District_Mapping_PRO.aprx")
m = aprx.listMaps("name_of_map")
inspector_layer = m.listLayers('name_of_layer')

arcpy.CopyFeatures_management(inspector_layer, r"C:\\testfolder\\name_of_layer.shp")

I am getting the following error: 'list' object has no attribute 'listLayers' on the line that reads inspector_layer = m.listLayers('name_of_layer'), hence my question of how to do it correctly.


texts:
I've been using GDAL to work with geospatial PDFs to make geotifs and have been running into some issues, so I started working with QGIS. A fine product I might add.

When I select a directory containing PDFs, the browser panel shows 2 PDFs listed for each actual PDF. One looks like a checkerboard and the other like a jellybean.

Dragging the jellybean to the work area pops up a window listing some layers from the PDF. When compared to the layers in the Acrobat view of the PDF, some are missing.

Dragging the PDF jellybean looking icon to the work area freezes the program. I believe that it is rendering.

Is there a way to get all layers to show and deselect the ones that I wish to remove?


texts: Field Calculate Degree Minute Second in Different Format?

texts:
I have a field that is in degree minutes seconds, populated from ET Geowizards COGO Inverse tool. It populates the values as 90-12-28.15 where 90 is the degree, 12 is the minutes and 28.15 is the seconds.

I need to use this field to label on a map. The format of the label is calling for the degree minutes seconds characters (i.e. 90° 12' 28.15") so I want to add a new field and use VB/Python to field calculate these characters in and remove the dashes.


texts: Geoserver SLD: setting fill color from property does not work

texts:
I'm trying to style a vector (Point) Layer in Geoserver 2.8 by using a property (attribute field) of the underlying  data source.
Here's the snippet:

<PointSymbolizer>
        <Graphic>
          <Mark>
            <WellKnownName>circle</WellKnownName>
            <Fill>
              <CssParameter name="fill">
                <ogc:PropertyName>HEXCOLOR</ogc:PropertyName>
              </CssParameter>
              <CssParameter name="fill-opacity">
                1
              </CssParameter>
            </Fill>
          </Mark>
        </Graphic>
</PointSymbolizer>

Where HEXCOLOR is a string like #fff000.
However, the point is rendered with white filling.

I have so far checked

the data (Q-GIS is able to render the color definitions)
the syntax


texts: wps-grass-bridge misses essential GRASS tool

texts:
I set up a Web Processing Service with 52°North and GRASS (using the wps-grass-bridge).
I get a similar result like in the online example here: 52°North WPS Demo
I assume it would be similar with PyWPS.

This is very nice. Now I would like to use i.cluster which needs to have the bands grouped in subgroups. In GRASS you would run i.group before running i.cluster. Unfortunately i.group is not available in any of the WPS. How can this be missing?

Does anybody have an idea how to make i.group available with the wps-grass-bridge? Or is there an alternative to i.group for defining groups on a raster dataset?


texts: Where to find 'Identify feature' tool in QGIS?

texts: It used to be an arrow with i. Now I have different version, I have gone through all the buttons and panels but cannot find it!


texts: How to model a geodesic grid in PostGIS?

texts:
Abstract
I've decided to use PostGIS for a hobby game project. It requires a lot of units allocated on the surface of a planet. I'd like to model not only the units (simple points), but also their areas of effect (vision, firing range etc.) and the planet's terrain using spatial data and queries.
The problem
I'd like to have a few sampled layers for things like the plant density. I've done a few GIS exercises on my university and observed that raster representation was most commonly used for such data.
The proposed solution and question
The problem with regular raster data is that it's best suited for smaller areas, which are almost flat, whereas I need a full sphere. I've learned about Geodesic Grid technique, which seems to fit my needs perfectly. However, I am not sure how to approach its creation and use in PostGIS. I've found this question, but it doesn't answer all of my needs.
In short, what I'd like to do is to create a model (with a chosen density) and define operations that would allow me to generate (fill in) the data in it, and then sample it back with some filtering.
Remark #1
I can imagine that some sort of clustering (variable density) models could be much more efficient in my scenario (some regions being much more important than the others, and thus needing more data). Would creating such a thing be much more difficult?
Remark #2
The Wikipedia article outlines the icosahedron tesselation procedure. I think I could implement it, but I am wondering if there are any caveats related.


texts: PyQGIS make screenshot of mapCanvas after setExtent is called

texts:
I have the following function:

def createScreenshots(layer):
  counter = 0
  for feature in layer.getFeatures():
    counter += 1
    point = feature.geometry().asPoint()
    qgis.utils.iface.mapCanvas().setCenter(point)
    qgis.utils.iface.mapCanvas().refreshAllLayers()
    qgis.utils.iface.mapCanvas().saveAsImage("D:\\\\m\\testing\\" + str(counter) + ".png")

Screenshots are not in real coordinates, .pgw header is not in right coordinates, therefore resulting png images are shifted.

Question: How to wait for rendering of canvas after setCenter is called?


texts:
I am writing a standalone application on QGIS3.8.  The environment is as follows:

@echo off
SET OSGEO4W_ROOT=C:\OSGeo4W64
SET QGIS=%OSGEO4W_ROOT%\apps\qgis
SET QGIS_PREFIX_PATH=%OSGEO4W_ROOT%\apps\qgis
CALL %OSGEO4W_ROOT%\bin\o4w_env.bat
CALL %OSGEO4W_ROOT%\bin\qt5_env.bat
CALL %OSGEO4W_ROOT%\bin\py3_env.bat
CALL %OSGEO4W_ROOT%\apps\grass\grass76\etc\env.bat
PATH %OSGEO4W_ROOT%\apps\qgis\bin;%PATH%
SET PYTHONPATH=%PYTHONPATH%;%OSGEO4W_ROOT%\apps\qgis\python
SET PYTHONPATH=%PYTHONPATH%;%OSGEO4W_ROOT%\apps\Python37
SET PYTHONPATH=%PYTHONPATH%;%OSGEO4W_ROOT%\apps\Python37\DLLs
SET PYTHONPATH=%PYTHONPATH%;%OSGEO4W_ROOT%\apps\Python37\Lib
SET PYTHONPATH=%PYTHONPATH%;%OSGEO4W_ROOT%\apps\Python37\Lib\site-packages
START "Visual Studio Code for QGIS" /B "C:\Program Files\Microsoft VS Code\Code.exe" %*

The code runs smoothly without any error. However, once I enter the debugging mode, the following error occurs. It looks like there is nothing wrong with my code, but something wrong with the debugging configuration. The tracback:

PS \Users\Richard\PycharmProjects\GPS_TT_Matrix_py3> cd 'c:\Users\Richard\PycharmProjects\GPS_TT_Matrix_py3'; ${env:PYTHONIOENCODING}='UTF-8'; ${env:PYTHONUNBUFFERED}='1'; & 'C:\OSGeo4W64\apps\Python37\python.exe' 'c:\Users\Richard\.vscode\extensions\ms-python.python-2019.6.22090\pythonFiles\ptvsd_launcher.py' '--default' '--client' '--host' 'localhost' '--port' '53041' 'c:\Users\Richard\PycharmProjects\GPS_TT_Matrix_py3\MainWin.py'
Traceback (most recent call last):
  File "c:\Users\Richard\.vscode\extensions\ms-python.python-2019.6.22090\pythonFiles\ptvsd_launcher.py", line 21, in <module>
    import ptvsd
  File "c:\Users\Richard\.vscode\extensions\ms-python.python-2019.6.22090\pythonFiles\lib\python\ptvsd\__init__.py", line 13, in <module>
    from ._vendored import force_pydevd
  File "c:\Users\Richard\.vscode\extensions\ms-python.python-2019.6.22090\pythonFiles\lib\python\ptvsd\_vendored\force_pydevd.py", line 37, in <module>
    'pydevd',
  File "c:\Users\Richard\.vscode\extensions\ms-python.python-2019.6.22090\pythonFiles\lib\python\ptvsd\_vendored\__init__.py", line 123, in preimport
    import_module(name)
  File "C:\OSGeo4W64\apps\Python37\Lib\importlib\__init__.py", line 127, in import_module
    return _bootstrap._gcd_import(name[level:], package, level)
  File "c:\Users\Richard\.vscode\extensions\ms-python.python-2019.6.22090\pythonFiles\lib\python\ptvsd\_vendored\pydevd\pydevd.py", line 19, in <module>
    from _pydev_bundle import pydev_imports, pydev_log
  File "c:\Users\Richard\.vscode\extensions\ms-python.python-2019.6.22090\pythonFiles\lib\python\ptvsd\_vendored\pydevd\_pydev_bundle\pydev_imports.py", line 8, in <module>
    import xmlrpclib
  File "C:\OSGeo4W64\apps\Python27\lib\xmlrpclib.py", line 194
    MAXINT =  2L**31-1
               ^
SyntaxError: invalid syntax

Any idea about how to fix this please?


texts: SQL script for postgis column name with colon

texts:
I have an sql script that creates some tables from the default planet_osm_point, planet_osm_line, planet_osm_polygon tables which are created with the osm2pgsql command line tool.

 DROP TABLE IF EXISTS "aero-poly";

 CREATE TABLE "aero-poly" AS (

  SELECT way,aeroway

  FROM planet_osm_polygon

 WHERE "aeroway" IN ('apron','runway','taxiway','helipad')

 ORDER BY z_order ASC );

CREATE INDEX "aero-poly_way_idx" ON "aero-poly" USING gist (way);

DROP TABLE IF EXISTS "agriculture";

CREATE TABLE "agriculture" AS (

 SELECT way,landuse FROM planet_osm_polygon

 WHERE "landuse" IN    ('allotments','farm','farmland','farmyard','orchard','vineyard')

ORDER BY z_order ASC );

 CREATE INDEX "agriculture_way_idx" ON "agriculture" USING gist (way);

and more ....

Some columns in the default tables mentioned above, have names that contain a colon (:) like this: seamark:type.

How may I use this name (seamark:type) in the sql script?

Is there a way to escape the colon (:)?


texts: Creating new layer in PyQGIS with same attributes of another preexisting and copy some elements in it

texts:
I'm new in QGIS and I'm trying to create my first plugin. I have a layer containing lots of elements of type Point, Line and Polygon. Suppose I want to create two new layers with the polygons above and below 15 vertices. I want the new two layers have the same attributes as the original one (CRS, colours, etc.) The general structure can be seen as

#original layer
layer = self.iface.activeLayer()
#HERE I WANT TO CREATE TWO EMPTY LAYERS WITH THE SAME ATTRIBUTES AS layer
layerAbove15 = NEW_LAYER(layerAbove15Name,same attributes as layer)
layerBelow15 = NEW_LAYER(layerBelow15Name,same attributes as layer)
#get layer features
feat = layer.getFeatures()
#loop over all elements
for elem in feat:
    #extract geometry
    geom = elem.geometry()
    #check if the element is a polygon
    if geom.type()==QGis.Polygon:
        #element as polygon
        poly = geom.asPolygon()
        #check the number of vertices of the outer ring
        if len(poly[0])>=15:
            #HERE I WANT TO COPY elem TO THE layerAbove15 LAYER
            layerAbove15.ADD(elem)
        else:
            #HERE I WANT TO COPY elem TO THE layerBelow15 LAYER
            layerBelow15.ADD(elem)

Is my pseudocode correct?


texts: How to install QGIS Server Ubuntu 12.04?

texts:
I am new user to Ubuntu (12.04) and am trying to set up QGIS server. I have QGIS desktop installed and working.

In regards to QGIS server I have used

apt-get install qgis-mapserver libapache2-mod-fcgid

and was successful but when I try the GetCapabilities request it fails and also when I try start/stopping apache2 I get the error "No Apache MPM package installed". Do I need to configure the apache2.conf file for this to work? Thanks for your time.

Below is example of error.

ben@ubuntu:~$ aptitude search apache2-mpm-
p   apache2-mpm-event               - Apache HTTP Server - event driven model
p   apache2-mpm-itk                 - multiuser MPM for Apache 2.2
p   apache2-mpm-prefork             - Apache HTTP Server - traditional non-threa
p   apache2-mpm-worker              - Apache HTTP Server - high speed threaded m
ben@ubuntu:~$ /etc/init.d/apache2 start
No apache MPM package installed


texts:
I am trying to update my feature class with some new parameters using IEditWorkSpace,
Following is my code ---

 IWorkspaceFactory workspaceFactory = new AccessWorkspaceFactory();
 IWorkspace workspace = new Workspace(workspaceFactory.openFromFile(inGDB, 0));

   IFeatureWorkspace ifw = (IFeatureWorkspace) workspace;
   IWorkspaceEdit iwe = (IWorkspaceEdit)workspace;
    iwe.startEditing(true);
    iwe.startEditOperation();

       IFeatureClass ifetclass =  ifw.openFeatureClass("CPWire");
       IFeatureCursor ifc = ifetclass.search(null, true);
       int field = ifc.findField("CREATIONUSER");
       System.out.println("CREATIONUSER :"+field);
       IFeature iff = ifc.nextFeature();
            while((iff = ifc.nextFeature()) != null)
            {
           iff.setValue(field, "xyz");
           iff.store();
       }
   iwe.stopEditOperation();
   iwe.stopEditing(true);
  }

While executing my code I get following error

Objects in this class cannot be updated outside an edit session

Why am I getting this error even if I have started the edit operation on workspace?


texts: Watershed Delineation for Specified Streams

texts:
I have a set of streams (currently a shapefile (line)) which I need to delineate the watersheds for. I have attempted to use various QGIS tools, however, each tool I have tried is for delineation of an entire region (i.e. DEM) rather than a few streams contained in a vector file.

My question is this: does anyone know of a tool designed to use vector files as an input (along with a filled dem) to generate watershed basin areas only for the lines in the file?

If more information is needed, please let me know and I will provide it and/or edit this post accordingly.


texts: How can I get the values returned by an OpenLayers bounds object in Longitude / Latitude format?

texts:
I've created a map and added a layer using the OpenLayers library:

map = new OpenLayers.Map('map');
var gphy = new OpenLayers.Layer.Google("Google Physical", {type: G_PHYSICAL_MAP});
map.addLayer(gphy);

The map loads and I can center my map:

map.setCenter(new OpenLayers.LonLat(-114.922485, 56.400224), 4);

I've then created a box handler, but the bounds returned aren't in Longitude / Latitude format. When drawing a box around the center of my map, my bound object returns the following:

left-bottom=(278,245) right-top=(304,211)

How can I convert these values into Longitude / Latitude?

The resolution thanks to unicoletti:

map.getLonLatFromPixel(new OpenLayers.Pixel(bounds.left, bounds.bottom))


texts: QgsVectorFileWriter.writeAsVectorFormat differents QGIS 2 QGIS 3

texts:
A long time ago I wrote a plugin for my own use. The plugin works fine in QGIS 2.18 but not in QGIS 3.1. I used the script 2to3 and solved most problems.The only problem I am left with is the following.

First my code that I use.

import os.path

CSVDIR = 'd:/QGIS/Temp_map/TestCSV/2018/'
DESTDIR = 'd:/QGIS/Temp_map/TestSHP/2018/'

if not os.path.exists(DESTDIR):
os.makedirs(DESTDIR)

for root, dirs, files in os.walk(CSVDIR):
    for file in files:
        if file.endswith('.csv'):
            fullname = os.path.join(root, file).replace('/', '/')
            filename = os.path.splitext(os.path.basename(fullname))[0]
            uri = 'file:///%s?crs=%s&delimiter=%s&xField=%s&yField=%s&decimal=%s&useHeader=no&' % (fullname, 'EPSG:4326', ';', 'Field_7', 'Field_6', ',')
            layer = QgsVectorLayer(uri, 'my_layer', 'delimitedtext')
            QgsVectorFileWriter.writeAsVectorFormat(layer, DESTDIR + filename + '.shp', 'CP1250', None, 'ESRI Shapefile')

    self.iface.messageBar().clearWidgets()
    self.iface.mapCanvas().refresh()

    QMessageBox.information(iface.mainWindow(), "CSV naar Shape",
        "Uw Shapefiles zijn aangemaakt!", QMessageBox.Ok )

This code works in QGIS 2.18 but give the following error in QGIS 3.0.

Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.0\apps\Python36\lib\code.py", line 91, in runcode
  exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 24, in <module>
  TypeError: QgsVectorFileWriter.writeAsVectorFormat(): arguments did not match any overloaded call:
  overload 1: argument 4 has unexpected type 'NoneType'
  overload 2: argument 4 has unexpected type 'NoneType'
  overload 3: argument 3 has unexpected type 'str'

What can I do to make this working?


texts: How to create a custom draggable rectangle overlay in Google Maps?

texts:
I'm trying to emulate a functionality in a certain app wherein a user can draw a rectangle and and hold & drag it. Then upon dragging, it shows the current length, width and area.

It's hard to put it in words so here's a GIF of what I'm trying to do.

Can this be achievable with google maps api alone or are there any third-party libraries for this functionality?


texts: Load LINESTRING in WKT-format from CSV to pgAdmin?

texts:
I've got a CSV file with a column containing LINESTRINGs that looks something like this:

1,"LINESTRING(20.331 54.33, 21.67 44.89, 12.67 13.5679)",2014-03-09 13:30:33,2014-03-09 13:31:33
2,"LINESTRING(25.331 14.33, 41.67 44.119, 15.67 23.5679, 34.55 48.11)",2014-03-10 13:30:33,2014-03-10 13:31:33
1,"LINESTRING(21.331 51.33, 11.67 44.32, 32.67 13.59)",2014-03-22 13:30:33,2014-03-23 13:31:33
5,"LINESTRING(10.331 24.22, 23.37 34.89, 22.17 12.29)",2014-03-24 13:30:33,2014-03-24 13:31:33

I created a table in pgAdmin 4 looking like this:

CREATE TABLE linestringtest (participant integer, geom geometry, startt time without time zone, endt time without time zone);

I now wanted to copy/import the CSV data, using the pgAdmin GUI, into the table but it doesn't work. I tried everything. Do I have to convert my linestrings into binary system first and if so how does it work?


texts: Extract values to point - values not being extracted

texts: I have a single band raster. When I use the information select tool I can see pixel values. I also have a shape file made up of points sitting over the raster. I want to extract the pixel values from the raster and have those values assigned to each point. I am using Extract values to Point in Arcmap 10.1.  The tool says it runs correctly and the output has a RASTERVALU column added. However all the values are zero, when they should be in the 800s.


texts: ArcGIS Pro Classify Raster Failure ERROR 002834

texts:
I am trying to run the Classify Raster tool in ArcGIS Pro 2.5 using an .ecd file generated from the ISO Classifier tool. I am using a composite of Sentinel bands run through the PCA tool and an additional SRTM layer to do the classification for the study area. The PCA layers, Composite image, and ISO Classifier ECD file generated just fine, but when I try to run the Classify Raster tool, I keep getting the 002834 error. I have tried to change the parameters, use different imagery as inputs into the tool, use alternative ECD files, but nothing seems to work. When I try to look up the error on the ArcGIS Pro website, there seems to be no documentation related to it.
How can I work around this error?


texts: Using euclidian distances as spacing for regular points grid in QGIS?

texts:
As a means of checking spatial autocorrelation in the selection of training samples for image classification, I am trying out semivariance analysis. I have managed to generate variograms from image pixels using gstat package in R. And I understand from this question (https://stackoverflow.com/questions/10769884/what-are-the-units-of-distance-in-gstat-variogram) that in my case the computed distance for the variogram is Euclidean distance. Now I wish to implement the range from the fitted variogram as spacing for regular points in QGIS, but I’m not sure how to translate Euclidean distance into point spacing.

Any Ideas, suggestions or experiences with what I want to do?
I have seen some publications that used this approach, with the distance in meters or pixel count but don't know how. I would also be open to try this in R or any other feasible software if pointed in the right direction..


texts: Wind Data for Newfoundland or New Brunswick

texts:
I am looking for wind data at an elevation of 100 metres for Newfoundland and New Brunswick that I can potentially bring into ArcGIS 10.2.2 to augment my wind turbine siting process.

Any information and links would be quite helpful.


texts: Convert raster values to categories using ArcGIS

texts:
I am currently processing the wetlands GIS data for North Carolina (obtained from https://www.fws.gov/wetlands/Data/Data-Download.html) using ArcGIS Pro. I downloaded this data as a shape file and converted it to a raster, in which each pixel represents a category e.g. lake, freshwater wetland, riverine, estuarine/marine wetland, etc. The issue is that when I clip this raster to any extent it replaces the categories with corresponding numerical values.

Is there a way to keep the categories or is there a way to relabel all the values as strings?

I attached an image of my raster symbology before and after clipping so you get a better idea of what my issue is:


texts: How to create a model which executes a SQL stored procedure with user-provided parameters?

texts:
How can I create an ArcToolbox that executes a SQL stored procedure and asks the user for two inputs?

As it stands now, I have to log in to my SQL server and run an execute statement for the stored procedure that requires two inputs from the user.

I'm trying to streamline my process to give people with no SQL experience the ability to generate these data through an ArcToolbox.


texts:
I would like to use Goode's Homolosine projection with GeoTools. This projection is fully supported by both PROJ and GDAL. However, when I try to create a CoordinateReferenceSystem object with this projection, GeoTools complains it has no transformation. For instance, the following little programme:

import org.geotools.referencing.CRS;
import org.geotools.referencing.crs.DefaultGeographicCRS;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;

public class WGS84toHomolosine
{
    static CoordinateReferenceSystem WORLDCRS = DefaultGeographicCRS.WGS84;
    static String wkt = "PROJCS[\"Homolosine\",\n" +
                        "    GEOGCS[\"WGS 84\",\n" +
                        "        DATUM[\"WGS_1984\",\n" +
                        "            SPHEROID[\"WGS 84\",6378137,298.257223563,\n" +
                        "                AUTHORITY[\"EPSG\",\"7030\"]],\n" +
                        "   AUTHORITY[\"EPSG\",\"6326\"]],\n" +
                        "        PRIMEM[\"Greenwich\",0,\n" +
                        "            AUTHORITY[\"EPSG\",\"8901\"]],\n" +
                        "        UNIT[\"degree\",0.0174532925199433,\n" +
                        "            AUTHORITY[\"EPSG\",\"9122\"]],\n" +
                        "        AUTHORITY[\"EPSG\",\"4326\"]],\n" +
                        "    PROJECTION[\"Interrupted_Goode_Homolosine\"],\n" +
                        "    UNIT[\"Meter\",1]]";
    static CoordinateReferenceSystem crs;
    static MathTransform transform;

    public static void main(String[] args)
    {
        // TODO Auto-generated method stub
        try
        {
            crs = CRS.parseWKT(wkt);
            transform = CRS.findMathTransform(WORLDCRS, crs, true);
        }
        catch (Exception e)
        {
            // TODO Auto-generated catch block
            e.printStackTrace();
        }

        System.out.println(transform.toWKT());
        System.out.println(transform.getClass().getCanonicalName());
    }
}

Regurgitates this output:

org.opengis.referencing.NoSuchIdentifierException: No transform for classification "Interrupted_Goode_Homolosine".
    at org.geotools.referencing.operation.DefaultMathTransformFactory.lambda$getProvider$1(DefaultMathTransformFactory.java:267)
    at java.util.Optional.orElseThrow(Optional.java:290)
    at org.geotools.referencing.operation.DefaultMathTransformFactory.getProvider(DefaultMathTransformFactory.java:264)
    at org.geotools.referencing.operation.DefaultMathTransformFactory.getDefaultParameters(DefaultMathTransformFactory.java:297)
    at org.geotools.referencing.wkt.Parser.parseProjection(Parser.java:632)
    at org.geotools.referencing.wkt.Parser.parseProjCS(Parser.java:988)
    at org.geotools.referencing.wkt.Parser.parseCoordinateReferenceSystem(Parser.java:232)
    at org.geotools.referencing.wkt.Parser.parseCoordinateReferenceSystem(Parser.java:212)
    at org.geotools.referencing.factory.ReferencingObjectFactory.createFromWKT(ReferencingObjectFactory.java:1068)
    at org.geotools.referencing.CRS.parseWKT(CRS.java:560)
    at org.isric.WGS84toHomolosine.main(WGS84toHomolosine.java:32)
Exception in thread "main" java.lang.NullPointerException
    at org.isric.WGS84toHomolosine.main(WGS84toHomolosine.java:41)

What is amiss here? Could there be something wrong with my system set up that would prevent GeoTools from finding PROJ's database?


texts: pgRouting - Route between multiple destinations - in order

texts:
I want to calculate the route between multiple points, respecting the order of the list.

For example, if I have the following list:

Lat: -40, Long: -30
Lat: -41, Long: -31
Lat: -40.5, Long: -30.5

I want the route to go from point 1 to 2, and then from point 2 to 3.

I don't want to calculate the best path from 1 to 2 and 3 (Travelling Sales Person) because that would give me a path 1-3-2, as the 1st point is closer to the 3rd point than the 2nd one.
I don't want to calculate the best path from 1 to 2, and from 1 to 3 (Mutliple destination Shortest Path Dijkstra), because that would give me a path from 1 to 2 and from 1 to 3.

I want it to calculate the best path from 1 to 2, and then from 2 to 3. All with just 1 query, using pgrouting.

It would be something like a TSRP but preserving the order of the points. Google already does that (https://developers.google.com/maps/documentation/javascript/directions#Waypoints) and if you want, they can optimize the order of the waypoints (using TSRP)


texts: Help dealing with a GeoTiff-file in ArcGIS

texts:
I'm doing data-analysis on satellite data and my source data is in GeoTiff-file format. What I would like to do with the data is:

1) Open the raster data in ArcMap etc.
2) Change the coordinate system from WGS84 --> ETRS-TM35FIN
3) Extract an area from the original data specified by ETRS-TM35FIN coordinates
4) Save the extracted area as plain .tif-file so that I can do data-analysis on it in Matlab

Can someone give the steps I need to do, in order to achieve my goals? I have no experience in ArcGIS software and the terminology etc. is very new to me, but I don't want to start reading a 500 page book just to get my data into a format so that I can analyze it in Matlab =)

The following picture will summarize my question:

Thank you for any help! =)

UPDATE:

I'm able now to extract specific area from the GeoTiff, but I can't save it as a plain .tiff-file:


texts: Welding line vertices or resume creating new line from end of old one in QGIS?

texts:
I'm creating a long and detailed (historical) border, and I guess it will take days.

How can I weld lines, or just resume creating a new one directly from the end of old fragment using QGIS?

When I click 'edit' and that small '+' at the end it just adds one(why?) segment and then stops without me right clicking it.

Edit: I've found 'Advanced Line Editor' plugin, but still it doesn't connect lines/vertices, however splitting works.


texts: Using ArcGIS conditional raster calculator?

texts:
I have a raster called "aspect" and I want to add 360 to every cells which have a current value lower than 100.

How can I do this using raster calculator?


texts: How can I display a coordinate grid in the map canvas?

texts: I want a live renderizing grid on the map canvas in QGIS, during processing, not on the map composer. Does anyone know how to do that?


texts: Find out ESRI REST URL for web map which is hiding it

texts:
I would like to extract the pipeline geometries from https://www.arcgis.com/home/webmap/viewer.html?webmap=a00c3b5cee4e4fe0b238b5e05ed80204

It shows the existing gas & oil pipelines crossing through Pennsylvania and was created by FracTracker for http://www.roverpipelinefacts.com/about/route.html . Despite it's a NGO they don't allow you downloading their data.

If I should know the REST FeatureServer URL I could download it (EXAMPLE):

wget "https://services.arcgis.com/jDGuO8tYggdCCnUJ/arcgis/rest/services/AtlanticCoastPipeline/FeatureServer/2/query?where=1=1&outFields=*&outSR=4326&f=geojson"

How can I find out this URL for the web map from above?


texts:
From the Google Engine scripts,

Is it available any Segmentation tool like the ArcMap "Segmentation Mean Shift" from ArcMap? I´m looking for crop-fields limit detection.
I tried the SNIC one explained by Gorelick at the 2018 Google Engine Summit (video, slides), but it needs the "Crop Data Layer" as an input, and I am located at Argentina. Maybe a SNIC workaround?


texts: Assigning a weight to thousands of values in the attribute table to then make a heat map in QGIS 3.10

texts:
I am working with some data that has about 19,000 values. This data has been imported as a text delimited file layer, as it is GPS (long and lat) point value data, that has an associated metadata or attribute field "hours". I have plotted the 19,000 points, and am trying to make a heat map based on the number of hours spent at each point. So, when I try the heatmap function and weigh the heat map by "hours" I do not receive the results I want. I understand that I can assign a weight to each of the 19,000 points based on how many hours spent at each point. I.e. The point with the highest amount of hours would be assigned the highest weight.
I have tried the CASE function in the field calculator but since this only works for individual values (as all the numbers of hours are different) and I have 19,000 values, I would like to know if there is a function where I can assign a weight for all 19,000 points based on how many hours spent, and then I can weigh the heat map based on the new attribute (the weights of the hours) that I just created.
Is this possible?


texts:
When requesting featureInfo from MapServer using getGetFeatureInfoUrl() of OL3 lib

var req_url = the_layer.getSource().getGetFeatureInfoUrl(
  evt.coordinate,
  map.getView().getResolution(),
  map.getView().getProjection(), //4326
  {
     'INFO_FORMAT': 'application/vnd.ogc.gml',
  }
);

it returns response in which the coordinates X and Y are in coorect order but after parsing the response theu are flipped:

$.ajax(req_url, {
            async: false
          })
          .done(function(resp){
            //console.log(resp);
            WMSfeatures = parser.readFeatures(resp);
            console.log(WMSfeatures);
           });

Coordinates in response from mapserver:

[[223.745,250.638],[232.028,255.22],[229.376,264.126]] (in EPSG:4326)

Parsed response:

[[250.638,223.745],[255.22,232.028],[264.126,229.376]] (in EPSG:4326)

What is wrong here?


texts:
I'm new to GIS and I'd like to apply some machine learning techniques to decent-quality satellite data for making certain types of predictions. This would potentially be for commercial purposes.

I'm looking for a way to access satellite imagery of the United States. I will need to select hundreds of high-resolution images to serve as training/testing data.

I was hoping, being so new to GIS, that there was either a free or inexpensive collection of satellite images of the U.S. to download, but its starting to sound like all the big companies (DigitalGlobe, etc.) that have to foot the bill to launch satellites into orbit aren't so keen on giving their data away to the everyman...

So I'm taking a wild swing here: are there any cheap-but-decent-quality data sets of satellite imagery, ideally of the U.S., that can be downloaded and even potentially used for commercial purposes?

Note: Google Earth/Engine is free to end users but if you read the fine print they don't want you using their images for commercial purposes.


texts: Convert multiple CSV files to shapefiles using Python?

texts:
At my internship, I have been asked to learn Python on the job.  It's been a rough road never having coded before so please forgive me if my problem is elementary.

The Lightning CSV files I am trying to convert have over 400,000 rows.
The headers go as follows:  mon, day, year, hr, min, sec, lat, lon, ht, type, Ip, Id

I found the answers to this question as potential solutions. I struggle to modify this code in order to convert the CSVs into shapefiles.

The path for the CSVs is C:\Users\zherran\Desktop\shp\enwest2015*.csv

The path where I want the shapefile C:\Users\zherran\Desktop\shp\shapefiles

import shapefile as shp
import csv

out_file = 'enwest3days.shp'

#Set up blank lists for data
mon,day,year,hr,min,sec,lat,lon,ht,type,Ip,Id=[],[],[],[],[],[],[],[],[],[],[],[]

# sample row mon,day,year,hr,min,sec,lat,lon,ht,type,Ip,Id
#            6,2,2015,7,27,16.6,41.5,-101.4,13789.5,1,-7900,26abb8a

#read data from csv file and store in lists
with open('enwest*.csv', 'rb') as csvfile:
    r = csv.reader(csvfile, delimiter=';')
    for i,row in enumerate(r):
        if i > 0: #skip header
            mon.append(int(row[0]))
            day.append(int(row[1]))
            year.append(int(row[2]))
            hr.append(int(row[3]))
            min.append(int(row[4]))
            sec.append(float(row[5]))
            lat.append(float(row[6]))
            lon.append(float(row[7]))
            ht.append(float(row[8]))
            type.append(int(row[9]))
            Ip.append(float(row[10]))
            Id.append(complex(row[11]))

#Set up shapefile writer and create empty fields
w = shp.Writer(shp.POINT)
w.autoBalance = 1 #ensures gemoetry and attributes match
w.field('X','F',10,8)
w.field('Y','F',10,8)
w.field('Date','D')
w.field('Target','C',50)
w.field('ID','N')
w.field('ID','N')
w.field('ID','N')
w.field('ID','N')
w.field('ID','N')
w.field('ID','N')
w.field('ID','N')

#loop through the data and write the shapefile
for j,k in enumerate(x):
    w.point(k,y[j]) #write the geometry
    w.record(k,y[j],date[j], target[j], id_no[j]) #write the attributes

#Save shapefile
w.save(out_file)


texts:
I am trying to adapt the openlayers code on using cql on the US states layer [topp:states] (http://dev.openlayers.org/examples/cql-format.html). I am able to change the wfs layer from the http://demo.boundlessgeo.com/geoserver/ version to the local http://localhost:8080/geoserver/ but the localised code result in an empty map (without the wfs ... only base layer shows up) if I change the following code (found in cql html right before the end of the html body):

src="http://demo.boundlessgeo.com/geoserver/wfs?service=WFS&version=1.0.0&request=GetFeature&typename=topp:states&outputFormat=text/javascript&format_options=callback:loadFeatures" type="text/javascript"

TO:

src="http://localhost:8080/geoserver/wfs?service=WFS&version=1.0.0&request=GetFeature&typename=topp:states&outputFormat=text/javascript&format_options=callback:loadFeatures" type="text/javascript"


texts: How to find the area of raster values within polygon zones?

texts:
I'm using ArcGIS 9.3 and I just can't seem to get this right - I've tried tabulate areas, zonal stats, converting to polygons, etc. I'm missing something somewhere.

I have a USGS impervious surface raster data set that has been clipped to my analysis zones (watersheds). I reclassified the values (100) into 10 values and then added fields and calculated the areas in sq mi and sq km throughout my analysis extent of each value. But what I really want is those areas calculated per analysis zone - 15 in all. I don't have classes here per se, just different levels of imperviousness. I want to see how much of each value is in each zone. It seems the tabulate areas tool should work, but I don't understand the numbers I'm seeing in the table. I'm choosing analysis name as my zone field and the sq mi as my class field. Any help would be much appreciated to a relatively new GIS user.

Thanks


texts:
I already ask how to print a map without screwing the "sync" between my map and the Openlayers background. I have tried 2.6.0 and 2.8.1 versions getting only deep frustration.
I have used the scale suggested to have optimal results, getting a wonderful series of failures.
On the screen the map looks perfect, once printed on PDF it slides out of sync.
I want to think I am doing something wrong because is hard to believe after so many years we have still to deal with this bug.
I love QGIS and what it allows me to do, but not been able to print my maps really make me @#$%$#@@#$##$@#@@#$&*%^%.
Help please.


texts: Get FeatureClass name from Layer properties in ArcObjects C#

texts:
I am newbie in ArcObjects. I am trying to fetch the featureclass name from layer properties under Data Source.
I am able to get the red-marked text via ILayer, IFeatureLayer, IFeatureClass. Is there any way to get the yellow-marked text, as it remained constant if I add same layer with different Definition Query?


texts: How does US Census map the input to shapes?

texts:
I've loaded a bunch of the census data. We're looking at using it for both geocoding of addresses, and also for census fact lookups (by tract or block or whatever).

The issue that we've bumped into is in geo coordinate accuracy... depending on the service we use, we've seen the same text address return coordinates up to 10 miles apart. We have plans for addressing the quality of it all, but it raised an interesting question...

The question is, when the US census data is being collected, how is it associated to the shapes (block/tract/whatever)... do they have a GPS in-hand and log the coordinates onto the piece of paper, or do they just log the address... in the case of the latter, is it known what service or data source they use for geocoding the recorded address?


texts: Exporting Delorme Street Atlas draw files into KML

texts:
(Using Delorme Street Atlas Plus USA 2015)

Is it possible to export a draw file into a KML file?  How do I do so?


texts: Finding interior angles of a triangle using PyQGIS

texts:
I need interior angles of a triangle (QgsTriangle object). It has angles() method. There is an example in documentation. But sometimes the method doesn't give correct results I expect.
For example:
from math import degrees

triangle = QgsTriangle( QgsPointXY(0, 0), QgsPointXY(1, 2), QgsPointXY(2, 0) )
angles = [ round(degrees(a), 2) for a in triangle.angles() ]
print(angles)

# Out: [26.57, 36.87, 26.57]

If I change the order of the first two points I get the correct result. But it depends.
triangle = QgsTriangle( QgsPointXY(1, 2), QgsPointXY(0, 0), QgsPointXY(2, 0) )
angles = [ round(degrees(a), 2) for a in triangle.angles() ]
print(angles)

# Out: [53.13, 63.43, 63.43]

I checked the source code of angles() method. It uses QgsGeometryUtils.angleBetweenThreePoints method.
I don't determine the order of the points. It comes from QgsGeometryobject. So how can I be sure that the angles() method always gives the correct angle values?


texts: Converting cluster of Points to Polygons using QGIS

texts:
I have an input point layer that is grouped into clusters, so I want to have each cluster of points converted into polygon.
This question How to convert a set of points to a polygon based on a specific field in QGIS (this question was removed) seems to come close to what I wanted by it hasn't been answered.
So basically, I need an algorithm in QGIS that will convert each cluster of points (based on attribute field) as seen below into a polygons?

I have tried "Point2One" plugin and it gave me wrong result seen below...

A copy of my data is found here.


texts:
The Problem:

I'm pretty new to leaflet and I'm trying to add a outline of the states over my OSM layer via WMS. I've found several tutorials online, some directly from the leaflet website but haven't had any success.

My code:

var map;
var feature;

function load_map() {
    map = new L.Map('map', {zoomControl: true});

    var osmUrl = 'http://<my_server>/osm_tiles/{z}/{x}/{y}.png',
            osmAttribution = 'Map data &copy; 2012 OpenStreetMap contributors',
            osm = new L.TileLayer(osmUrl, {maxZoom: 18, attribution: osmAttribution});

    map.setView(new L.LatLng(0, 0), 3).addLayer(osm);

    var states = L.tileLayer.wms('http://<my_server>:8080/geoserver/States/wms', {
    format: 'img/png',
    transparent: true,
    layers: "null"
    });

    states.addTo(map);
}

window.onload = load_map;

When I open the console in google chrome there are no reported errors with my java script or connecting to the server. When the map is displayed it has no problem displaying the OSM layer, but the WMS layer is missing.


texts:
I just install "SFCGAL". When I use "ST_3DIntersection", It get below error:

ERROR: attempt to redefine parameter "postgis.backend" SQL state: XX000

How can I fix?

I already saw http://postgis.net/docs/postgis_installation.html, But I do not know how should I do.


texts: Determining if shapefile has elevation or z coordinate values stored using ArcGIS Desktop?

texts:
I have a water main (line feature) that I wanted to get the elevation. As this is a subsurface feature I can't use any surface file to find these calculations. I wanted to use the calculate geometry tool to find the Z coordinate but that is not available.

It will allow me to calculate X & Y for start, mid pt, and end pt, so I'm assuming that this shapefile does not have elevation information stored in it?

How can I confirm that this shapefile only has X-Y information stored?


texts:
I'm wondering about the best strategy to incorporate user comments and user added points to a map. The best example that I can find of this is ushahidi with an example of how users can add information here:

I understand how to create this through the use of a web-form and a WFS-T type approach so that a user can add a point which can be stored in a postgresql/postgis database. Where I would appreciate advice is on a good strategy for displaying user-generated data back to the users, so that additional comments and interactions can be achieved.

I've looked through the source code of ushahidi and find it more sophisticated than I can handle at the moment. Are there any other examples or approaches of collaborative web-maps out there that illustrate good ways of doing this?


texts: Iterating through an numpy array and index to a value in another numpy array

texts:
I am struggling to get this code to work I want to iterate through an numpy array and based on the result, index to a value in another numpy array and then save that in a new position based on that value.

     # Convert the sediment transport and the flow direction rasters into Numpy arrays
    sediment_transport_np = arcpy.RasterToNumPyArray(sediment_transport_convert, '#', '#', '#', -9999)
    flow_direction_np = arcpy.RasterToNumPyArray(flow_direction_convert, '#', '#', '#', -9999)

    [rows,cols]= sediment_transport_np.shape
    elevation_change = np.zeros((rows,cols), np.float)

   # Main body for calculating elevation change

    # Attempt 1
    for [i, j], flow in np.ndenumerate(flow_direction_np):
        if flow == 32:
            elevation_change[i, j]  = sediment_transport_np[i - 1, j - 1]
        elif flow == 16:
            elevation_change[i, j] = sediment_transport_np[i, j - 1]
        elif flow == 8:
            elevation_change[i, j] = sediment_transport_np[i + 1, j - 1]
        elif flow == 4:
            elevation_change[i, j] = sediment_transport_np[i + 1, j]
        elif flow == 64:
            elevation_change[i, j] = sediment_transport_np[i - 1, j]
        elif flow == 128:
            elevation_change[i, j] = sediment_transport_np[i - 1, j + 1]
        elif flow == 1:
            elevation_change[i, j] = sediment_transport_np[i, j + 1]
        elif flow == 2:
            elevation_change[i, j] = sediment_transport_np[i + 1, j + 1]

    #sediment_flow_np = np.concatenate((flow_direction_np, sediment_transport_np), axis = 0)
    #case1 = numpy.where((flow_direction_np == 32)) & (sediment_transport_np > 0),

    # Convert the Numpy arrays back to rasters
    #sediment_transport_raster = arcpy.NumPyArrayToRaster(sediment_transport_np, bottom_left_corner, raster_cell_width, raster_cell_height, -9999)
    #flow_direction_raster = arcpy.NumPyArrayToRaster(flow_direction_np, bottom_left_corner, raster_cell_width, raster_cell_height, -9999)
    elevation_change_raster = arcpy.NumPyArrayToRaster(elevation_change, bottom_left_corner, raster_cell_width, raster_cell_height, -9999)
    elevation_change_raster.save(output_raster)

The error i get is:

Running script elevation_change...

Traceback (most recent call last):
  File "", line 606, in execute
IndexError: index (655) out of range (0<=index<655) in dimension 0

Failed to execute (elevation_change)


texts: Programmatically setting key value pairs for ValueMap widget of QGIS?

texts:
Using QGIS 2.6 and Python 2.7 on Windows 7, I am trying to programmatically insert key value pairs into a ValueMap widget via a csv file. Issue is that some of the values are properly inserted and formatted, while some are not.

This csv value is formatted like so:

            ;
20          ;              20
25          ;              25
32          ;              32
40          ;              40
50          ;              50
60          ;              60
63          ;              63
75          ;              75
80          ;              80
90          ;              90
100         ;              100
110         ;              110
125         ;              125
140         ;              140
150         ;              150
160         ;              160
200         ;              200
225         ;              225
250         ;              250
300         ;              300
315         ;              315
350         ;              350
400         ;              400
500         ;              500
1/4''       ;              1/4''
3/8''       ;              3/8''
1/2''       ;              1/2''
3/4''       ;              3/4''
1''         ;              1''
1'' 1/4     ;              1'' 1/4
1'' 1/2     ;              1'' 1/2
2''         ;              2''
2'' 1/4     ;              2'' 1/4
2'' 1/2     ;              2'' 1/2
3''         ;              3''
3'' 1/2     ;              3'' 1/2
4''         ;              4''
5''         ;              5''
6''         ;              6''

I then created a dictionary from these values, programmatically set the widget, and configured the widget using this code:

resLayer = qgis.utils.iface.activeLayer()
with open('mycsv.csv') as f:
    f.readline()
    diamResDict = dict(csv.reader(f, delimiter=';'))
diamidx = resLayer.fieldNameIndex('Diametre')
resLayer.setEditorWidgetV2(diamidx, 'ValueMap')
resLayer.setEditorWidgetV2Config(diamidx, diamResDict)

The dictionary prints this:

{'140\t': '\t140', '250\t': '\t250', '50\t': '\t50', '\t': '\t\t\t\t', "1'' 1/4\t": "\t1'' 1/4", '150\t': '\t150', '40\t': '\t40', '125\t': '\t125', '400\t': '\t400', "3/4''\t": "\t3/4''", "1'' 1/2\t": "\t1'' 1/2", "3'' 1/2\t": "\t3'' 1/2", '25\t': '\t25', '100\t': '\t100', "1/2''\t": "\t1/2''", '300\t': '\t300', "4''\t": "\t4''", "6''\t": "\t6''", '60\t': '\t60', "5''\t": "\t5''", '80\t': '\t80', '32\t': '\t32', '225\t': '\t225', '200\t': '\t200', '110\t': '\t110', '160\t': '\t160', '350\t': '\t350', "1/4''\t": "\t1/4''", "3''\t": "\t3''", '90\t': '\t90', '500\t': '\t500', "2''\t": "\t2''", "1''\t": "\t1''", '315\t': '\t315', "3/8''\t": "\t3/8''", "2'' 1/2\t": "\t2'' 1/2", '20\t': '\t20', '63\t': '\t63', '75\t': '\t75', "2'' 1/4\t": "\t2'' 1/4"}

The result is this:

It's the apostrophes instead of keys that is causing me trouble.


texts: Leaflet OSM Cross origin issue when zoom level > 19

texts:
I am playing around a bit with Leaflet and run into an issue.

I am using the OSM tile service https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png.

But I get CORS problem when I zoom in closer than level 19. Why do I get it only if I zoom closer than 19? Up to 18, it works fine. Shouldn't I get CORS issues on all zoom layers if it's truly a CORS problem?

Full error message

  Access to Image at
  'https://a.tile.openstreetmap.org/20/554563/283403.png' from origin
  'http://localhost:4321' has been blocked by CORS policy: No
  'Access-Control-Allow-Origin' header is present on the requested
  resource. Origin 'http://localhost:4321' is therefore not allowed
  access. The response had HTTP status code 404.


texts: am interested in mapping tree canopy cover per unit area basically looking at how density has changed over time using QGIS. some one with answeres kindly share


texts: How to add my own GRASS algorithms to Processing?

texts:
I need to use algorithm r.bioclimatic from GRASS in QGIS but I see that this algorithm don't exist in Processing Toolbox. A search a way of add this algorithm to Processing and in GRASS and Sextante page if find this way:

HOWTO adding own GRASS algorithms to Sextante (simply by adding a new description file)

I create a description file according to the rules but now I don't know how to run and if this will run in GRASS.

The first parameters of r.bioclimatic (tmin, tmax....) are columns of tables? I don't know which type of input I have.


texts: Software Recommendations: make a flyover video from GPX log?

texts: I m trying to make a video similar to this (at 53 seconds): https://youtu.be/OEU5BCZV6Sc?t=53 from GPX logs we have from our drone. Google Earth Pro software doesn't allow nearly enough customization, but I'm new to all of this and haven't been able to find what I should use. Paid software is fine as well.


texts: Drawing concentric circles with different radii in QGIS

texts: I have table that contain point data and multiple distance value for the point data in three different direction. I would like draw arc of 65degree in each of these three directions where the radius of the arc will vary according to a given distance data. It could  be in one direction I may have 5 arcs as only 5 distance values are non-zero and in another direcion i will have only 3 arcs as 3 values are non-zero. Is this possible ?


texts: Mosaic images in double precision pixel type

texts:
I am trying to mosaic the raster images, but those images are in the double-precision pixel type and their pixel depth is 64 bit.
How can I mosaic that raster image?
I am using ArcMap 10.2 version.
Actually, I need to mosaic the image with the 32-bit float but I have an image with pixel type double precision and I don't have knowledge about that type.


texts: Using field calculator if statement to update another field using another field's values?

texts:
Using VBA in field calculator I am trying to figure out the expression that if 3 conditions are met populate a field using another field's value and I cant get it to work! CALC_CLASS is a joined table and SPEED LIMIT is a Double type, MMS_TEMP is text.

I am using VBA and ARC 10.3

if [SPEED_LIMIT]=>91 and [SPEED_LIMIT] <= 100 and [MMS_TEMP]="53,000 or more" then
RD_BASELINE.CLASS_MMS= [CALC_CLASS_MMS.PSL_91_100]

elseif [SPEED_LIMIT]=>51 and [SPEED_LIMIT] <= 60 and [MMS_TEMP]="53,000 or more" then
RD_BASELINE.CLASS_MMS= [CALC_CLASS_MMS.PSL_51_60]

else
RD_BASELINE.CLASS_MMS=[CALC_CLASS_MMS.PSL_1_40]

end if**


texts: Raster area grouped by class and by a second raster

texts:
I have two categorical rasters, and I want to get the area of each category on the first further grouped by each category on the second in R, I first tried with zonal stats but you only get the (max, mean, min, count) of the values of the first raster, so now I'm trying to use freq and then multiply by the area the pixel resolution yields, is this ok? Us there a better way to achieve it?

library(raster)
r <- raster(ncols = 5, nrows = 6)
r[] <- sample(rep(1:5, 6), 30)
y <- raster(ncols = 5, nrows = 6)
y[] <- sample(rep(1:3, 10), 30)

df <- data.frame()
for(i in 1:4) {
  fkwence <- freq(r*(y == i))
  fkwence <- as.data.frame(fkwence)
  fkwence$y_raster_value <- i
  df <- rbind(df, fkwence)
}

and the outcome looks like this:

df
   value count y_raster_value
1      0    20              1
2      1     2              1
3      2     2              1
4      3     2              1


texts: Connect Multipoint Geometry type to create line

texts: As you can see in the pic below i have some points "Sông Vĩnh Trà". Their Geometry type are Multipoints. They are parts of a river, so i want to connect them to display the river on the map. I tried to search for "connecting Multipoint arcgis" but i still can't find solution.


texts: OpenLayers are not rendering MapServer layer

texts:
I have a MapServer layer, Roads in this instance, and I've set its STATUS to ON. I want to render this as a WMS service using OpenLayers. So I've set all the Metadata Section in the LAYER section in the MAP file. Then in OpenLayers I have the following snippet:

 var params = {
              LAYERS: 'Roads',
              FORMAT: 'image/png',
              VERSION: '1.1.1',
       }

       var source = new ol.source.ImageWMS({
            url: 'http://localhost/cgi-bin/mapserv?map=/home/user/mapdata/fifth.map&layer=Roads',
          params: params
       });

       var layer = new ol.layer.Image({
            source: source
       });

       var center = ol.proj.transform([21.768006, 39.556202], 'EPSG:4326', 'EPSG:3857');

       var view = new ol.View({
              center: center,
              zoom: 15
       });

       var map = new ol.Map({
            target: 'map',
            view: view,
            renderer: 'dom'
       });

       map.addLayer(layer);

However, the layer is never rendered, although in the url I pass the layer=Roads parameter.
If I set the STATUS to default in the map file everything works well but in this occasion, not. Any idea?


texts: Plotting large numbers of objects on a chart using R and WKT format

texts:
I am retrieving thousands of geometry lines in WKT format from a Microsoft SQL Server database.  I would like to plot all of those lines in one map, but I am having some difficulty with the code.  How should I modify the following code to reach my goal?

library(RODBC)

png(file="examplex.png", width=600, height=480)
setwd("C:/ArcR")

con <- odbcDriverConnect('driver={SQL Server};server=SERVERNAME;database=DBNAME;trusted_connection=true')

objects_1 <- sqlQuery(con, 'SELECT TOP (1) Shape.STAsText() as ShapeWKT FROM TABLENAME ;')

things <- vector("list", 1)

z = 0
for(line in objects_1$ShapeWKT)
{
  {
    things[[z+1]]<-readWKT(line)
  }
  z = z + 1
}

plot(things[[1]])

dev.off()

Yes, the code does refer to only one returned record, but I wanted to start the question with code in working order, for those who may not know the answer, but are still capable of testing possibilities.


texts: OpenLayers-5 reproject OSM in a custom projection

texts:
Following OpenLayers examples, I try to reproject OSM in a custom projection. The view is empty, no OSM tile is downloaded.

index.js :

import {get as getProjection} from 'ol/proj';
import Map from 'ol/Map';
import 'ol/ol.css';
import OSM from 'ol/source/OSM';
import proj4 from 'proj4';
import {register} from 'ol/proj/proj4';
import TileLayer from 'ol/layer/Tile';
import TileWMS from 'ol/source/TileWMS';
import View from 'ol/View';
import * as Gp from 'geoportal-extensions-openlayers';

proj4.defs('EPSG:2154','+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3' +
'+x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs');
register(proj4);
var franceProjection = getProjection('EPSG:2154');

 var couches = [
  new TileLayer({
  source: new OSM()
})
];
var map = new Map({
 target: 'map',
 layers: couches,
 view: new View({
projection: 'EPSG:2154',
center: [524326, 6700656],
zoom: 14
 })
});

package.json

{
"name": "iliad_ol5.3.0",
"version": "1.0.0",
"main": "index.js",
"scripts": {
"test": "echo \"Error: no test specified\" && exit 1",
"start": "parcel index.html",
"build": "parcel build --public-url . index.html"
},
"keywords": [],
"author": "",
"license": "ISC",
"dependencies": {
"geoportal-extensions-openlayers": "^3.0.7",
"ol": "^5.3.0",
"proj4": "^2.6.0"
},
"devDependencies": {
"parcel-bundler": "^1.12.4"
},
"description": ""
}

I know that it is a basic projection issue but even that's hard for me.


texts: Replicating ArcMap's Point Sketch Geometry Editing (Feature Templates) in ArcObjects

texts:
I'm looking for a way to place, select, move, delete features inside a custom C# ArcObjects application (using AxMapControl). I found the exact behavior in ArcMap, but I want to save the user some steps and avoid all of the dialogs, etc. It appears to be creating a point sketch geometry.

In ArcMap, I can do the following steps:

Open a new map doc
Choose start editing from the Editor toolbar
Then I get a dialog asking me to choose the layer to edit
After selecting a layer, a panel appears that is titled "Create Features". Under it, there is another panel named "Construction Tools"
When I click on a feature in the upper panel, I get a "Point" option in the lower panel
Now when I move the mouse over the map, the cursor has the same icon as the feature I selected.
Clicking will place a feature, then I can drag or delete.

So ... how can I replicate this behavior in a custom app using ArcObjects. I've already been able to replicate the behavior of other tools by doing something like this:

ICommandPoolEdit _commandPool = new CommandPoolClass();
ICommand solveRouteToolCmd = new ControlsNetworkAnalystSolveCommand();
_commandPool.AddCommand(solveRouteToolCmd, null);

ICommand command = commandPool.get_Command((int)solveRouteToolCmd );

//pass the MapControl to the oncreate method
command.OnCreate(_axMapControl.Object);
//run the tool
command.OnClick();

However, I can't seem to find a tool class object that corresponds to this behavior. In this example ControlsNetworkAnalystSolveCommand class is the hook to the same behavior that happens when you click the solve route button in ArcMap. I'm looking for that same hook for the scenario described above.

Note that I found the tools by looking at the customize menu dialog. They are in the commands tab, in the editor category. The specific command is called Point and the description says, "Creates a point sketch geometry".


texts: What unit of measurement is rxwaveform from the GEDI L1B dataset in?

texts:
Rxwaveform is listed in GEDI's data level product L1B's data dictionary as "The corrected receive (RX) waveforms. Use rx_sample_count and rx_sample_start_index to identify the location of each waveform.
" This can also be seen in the attached image. I understand that waveforms are typically measured by amplitude as a function of time(or height), but are the values listed in rxwaveforms amplitudes? It doesn't explicitly say amplitude in the data dictionary, and what unit of measurement is amplitude in (if it is indeed amplitude). The second image is the data values of rxwaveform in matlab; and what is the significance of rx_sample_count in relation to rxwaveform?


texts: QGIS color picker works only inside application

texts:
Tried to use QGIS Color picker for picking color for layer outside the QGIS window.
When Color picker is out of QGIS desktop application window, it changes to a simple cursor. Windows 8 x64, Windows 7 x86, x 64). However, everything works fine on my notebook (Windows 7 x64).
QGIS version - 2.18.


texts: I can apply tasseled Cap transformation in surface reflectance (not ToA)?

texts: I can not understand if I should adjust the Tasseled Cap index with the surface reflectance values or the Top of Atmosphere


texts: Extracting Mean from Raster and Writing to Table in ArcPy?

texts:
I am simply trying to extract the mean of each raster within a geodatabase and write the value to a table. I have done this successfully with another set of rasters, but now I keep getting the following error: Parsing error SyntaxError: invalid syntax (line 8). I am at a loss for what is causing this.

import arcpy
from arcpy import env
env.workspace = "G:/Research/Sand Dams/Datasets/Ukraine.gdb/"

rasters=arcpy.ListRasters("*", "All")

for inRaster in rasters:
    UkrMEANResult = arcpy.GetRasterProperties_management(inRaster, “MEAN”)
    UkrMean = UkrMEANResult.getOutput(0)
    arcpy.CalculateField_management("RasterList","UkraineMean",UkrMean)


texts: Converting coordinates (latitude,longitude) to X,Y using PostGIS?

texts: How can I convert the latitude and the longitude coordinates of a rectangle in Cartesian coordinates, so that I can store them in a geometry field using PostGIS?


texts: Using arcgisimage along with python basemap

texts:
I am utilizing maps provided through the ArcGIS rest service when i plot information in Python, but have some problems understanding the different options I have when I request I map. The documentation is really really sparse.

First, here is my code:

map = Basemap(epsg=3395,
llcrnrlon=27.692793, llcrnrlat=70.390829,
urcrnrlon=28.728254, urcrnrlat=70.855549)
map.arcgisimage(service='World_Topo_Map', xpixels=4000, xsize=16000, verbose= True)

This, for some reason leaves me with a topographic map where only the bottom half have good resolution, while the rest is really lacking. Anyone got experience in using this, and know what I do wrong?

This is the resulting http call (from the function above), and you can also see what I describe (half of the image is high resolution, half is low). Warning, large file.

http://server.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer/export?bbox=3082747.6154,11156828.0966,3198014.60666,11312643.9336&bboxSR=3395&imageSR=3395&size=8000,8000&dpi=96&format=png32&f=image