texts: Minimum Information to be kept on importing OSM data to use it for creating random start and ent points

texts:
I want to only import the minimum of data from OpenStreetMap, to be able to distribute random start and end points for a simulation in administrative areas in germany. I found this file https://github.com/openstreetmap/osm2pgsql/blob/master/empty.style which will act as my starting point.

Basically what I need are the bundaries of the districts, what type of distict it is and the streets. Nice to have would be if it is an area, where people life, office buildings or factories, since people commute between their home and work.

What I came up with so far is this. My Question to you: Am I missing something? Did I put in too much?

Generally I want to have a small database with the correct datatypes etc.

# osm2pgsql style file for simulation of commuter routes.
# The data is mainly used to determine random start and destination points
# for the route generation

# OsmType   Tag                     Type    Flags
# Insert your own columns here, or change phstore to polygon below
node,way    access                          text    linear
node,way    admin_level                     int4    linear
node,way    amenity                         text    polygon
node,way    area                            text    # hard coded support for area=1/yes => polygon is in osm2pgsql
way         boundary                        text    polygon
node,way    building                        text    polygon
way         de:amtlicher_gemeindeschluessel int4    polygon
way         de:regionalschluessel           int4    polygon
node,way    highway                         text    linear
node,way    landuse                         text    polygon
node,way    name                            text    linear
node,way    name:de                         text    phstore
node,way    office                          text    phstore
node,way    place                           text    polygon
node,way    population                      text    linear  # sadly text is used here as well
node,way    route                           text    linear
node,way    type                            text    linear
node,way    z_order                         int4    linear  # This is calculated during import
way         way_area                        real            # This is calculated during import

# Deleted tags
node,way    note                    text    delete
node,way    note:*                  text    delete
node,way    source                  text    delete
node,way    source_ref              text    delete
node,way    source:*                text    delete
node,way    attribution             text    delete
node,way    comment                 text    delete
node,way    fixme                   text    delete
node,way    created_by              text    delete
node,way    odbl                    text    delete
node,way    odbl:note               text    delete
node,way    SK53_bulk:load          text    delete
node,way    tiger:*                 text    delete
node,way    NHD:*                   text    delete
node,way    nhd:*                   text    delete
node,way    gnis:*                  text    delete
node,way    geobase:*               text    delete
node,way    accuracy:meters         text    delete
node,way    sub_sea:type            text    delete
node,way    waterway:type           text    delete
node,way    KSJ2:*                  text    delete
node,way    yh:*                    text    delete
node,way    osak:*                  text    delete
node,way    kms:*                   text    delete
node,way    ngbe:*                  text    delete
node,way    naptan:*                text    delete
node,way    CLC:*                   text    delete
node,way    3dshapes:ggmodelk       text    delete
node,way    AND_nosr_r              text    delete
node,way    import                  text    delete
node,way    it:fvg:*                text    delete

EDIT: Updated the style file


texts: How to add label in one direction on cross-sections?

texts:
I would like to add cross-section labels in one direction to the mouth of the river. Unfortunately, in some places label turns 180 degrees as pictured. Do you know some automatic way to solve this problem? I tried to use the advanced label settings (label manager) but it failed.

XS LABEL http://zibimobile.pl/adik/xs_label.jpg


texts: What is the "Auto Open Form" checkbox used for?

texts:
Windows 7
QGIS 2.6 - 64 bit
Feature - Identify Results

Could someone help me understand what is "Auto Open Form" checkbox used for.


texts: Selecting Specific Elevation Areas Contiguous to the Coastline for Coastal Flood Mapping

texts:
This is less a problem than an inquiry into a better way of doing something. I am curious how people might go about selecting, say, all areas less than or equal to 10 feet (from a DEM) that are contiguous to the coastline. The purpose of this exercise would be to map low-lying areas along the coast that may flood due to surge during a coastal storm.

Here is an abridged version of how I've been doing this sort of thing:

1) Extract from DEM all values less than or equal to 10 feet using the (raster) clip tool.
2) Convert the raster to a polygon using the Raster to Polygon tool.
3) Select areas contiguous to the coastline layer:
      a) Select all polygon areas that intersect the coast line. Export selection.
      b) Add to selection all areas that intersect the new selection. Export.
      c) Repeat this process until the difference in the number of attributes from one exported selection           to the next reaches 0, meaning, you are no longer selecting new features when you select            intersecting features.

4) Delete the remaining "orphaned" areas.

It's really step "3" that I'm trying to find an alternate method to. It can take up to 20 iterations of this step to select these areas for NYC. It just seems so tedious. I'm really curious how others might go about performing this task.

Update: I've come up with another method that eliminates the very iterative Step 3 from above. I am still curious to read how others might accomplish this task though.

Alternative Method:

1) Extract from DEM all values less than or equal to 10 feet using the (raster) clip tool.
2) Give every raster cell the same value (i.e. create a new field and have all values in that field = 1).
3) Convert raster to polygon using the Raster to Polygon tool.
4) Select and export features that intersect the coastline.
The above will get you your contiguous "flooded" area, but it will not have elevation values assigned to it, like the first step of steps above. To get those back, you'll need to do the following:

         A) Use the resulting features (shapefile or feature class) as a mask to clip the raster.
         B) Convert raster to polygon.


texts: Displaying Feature Layer Field Values on webpage without Map using ArcGIS API for JavaScript

texts:
I am just learning JavaScript for this mini-project.
I am looking to display the field values from a service feature layer on a web page using HTML and JS without loading a map and without any input or action from the user.
Service Feature Layer: https://services3.arcgis.com/GVgbJbqm8hXASVYi/arcgis/rest/services/Trailheads/FeatureServer/0
FieldID = TRL_NAME
What I've been working with so far is below from reading through tutorials and forums, but I can't get it working:
<!DOCTYPE html>

 <link rel="stylesheet" href="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.6/css/bootstrap.min.css" />

 Trails!

 Trail Name


texts: arcpy insert cursor problem - adding multiple polygons

texts:
I am reading rows from a database with a field "Coords" that contains polygon coordinates in the format 40.090170,-83.007704|40.089878,-83.007724..... and inserting each row as a polygon.

The first feature is drawn correctly, the second feature is drawn but  is connected to and duplicates the first feature, and the third contains all three and so on.

The coordinates are cumulative.  I tried deleting variables between iterations, but that didn't help.  What am I missing?

import pypyodbc
import arcpy
GdbFeat = r"\\server\Arc_Server_data\publish.gdb\OUPS"
spatialRef = arcpy.Describe(GdbFeat).spatialReference
conn = pypyodbc.win_connect_mdb(r"C:\GIS\db\COWOUPSx.mdb")
cursor = conn.cursor()
cursor.execute("SELECT TicketDate, TicketNum, Coords FROM Tickets WHERE (((TicketDate) = #11/15/2014#))")
rows = cursor.fetchall()
cursor.close
del cursor
conn.close
del conn
if rows:
        polygonArray = arcpy.Array()
        for row in rows:
                Coords = row[6]
                coordinatePairList = Coords.split("|")
                for coordinatePair in coordinatePairList:
                        coordinates = coordinatePair.split(",")
                        currentPoint = arcpy.Point(coordinates[1],coordinates[0])
                        polygonArray.add(currentPoint)

                polygon = arcpy.Polygon(polygonArray, spatialRef)
                polygonArray.removeAll
                with arcpy.da.InsertCursor(GdbFeat, ("Shape@", "TicketDate", "TicketNum")) as Acursor:
                        Acursor.insertRow((polygon, row[0], row[1]))

                del Coords
                del polygon
                del coordinatePairList
                del coordinatePair
                del coordinates
                del currentPoint

        del row
        del rows


texts: Converting tab to shapefiles with QuickImport gives fatal error?

texts:
I'm trying to convert tab to shape files.  When I run the code, I see the tab files(list) but I get the error using QuickImport.. Can somebody suggest a fix? My tab files are located in folders undernealth archive.

import arcpy
import env
import os

workspace = "Z:/STAGING/archive"

#returns @string list of feature classes
walk = arcpy.da.Walk(workspace, datatype="Any", type="All")

for dirpath,dirnames, filenames in walk:
    print filenames
#Set local variables for QuickImport
tabfiles = "Z:/STAGING/archive/**/*.TAB"
output_gdb = "c:/Temp/results.gdb"

# Execute Quick Import
arcpy.QuickImport_interop(tabfiles,output_gdb)

  Traceback (most recent call last):   File "C:\polygeo\Scripts\tab to
  shape.py", line 42, in
      arcpy.QuickImport_interop(tabfiles,output_gdb)   File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\interop.py", line 104, in
  QuickImport
      raise e ExecuteError: NULL pointer check failed.   File: geodbfeatureclass.cpp   Line #: 985 A fatal error has occurred. Check
  the logfile above for details Tool execution failed.

Failed to execute (QuickImport).

This is a minimal code snippet that I think will illustrate the same error:

import arcpy
tabfiles = "Z:/STAGING/archive/**/*.TAB"
output_gdb = "c:/Temp/results.gdb"
arcpy.QuickImport_interop(tabfiles,output_gdb)


texts:
Is it possible to calculate oblique photo scale by given vertical photo scale ?
I mean is there any kind formula relation between them or maybe is it possible to calculate rough approximation of oblique scale ?


texts: Merge raster assigning to one the overlapping area QGIS

texts:
I want to merge 2 raster files, but assigning the overlapped area to a specific one.

How could I do it?


texts:
What does the SyncOnce tag mean?

I have a Python toolbox and this xml document was created.

toolbox.tool.pyt.xml

<?xml version="1.0" encoding="UTF-8"?>
<metadata xml:lang="en">
    <Esri>
        <CreaDate>20180824</CreaDate>
        <CreaTime>13413800</CreaTime>
        <ArcGISFormat>1.0</ArcGISFormat>
        <SyncOnce>TRUE</SyncOnce>
    </Esri>
</metadata>

or is there a list of what all tags mean or information about the xml structure?

I know there are other tags that can go in this document like a summary tag which can be accessed/edited through the Item Description. However there is no edit section for SyncOnce.


texts: What are some good introductory books or articles about Open Source GIS for students coming from an ESRI background?

texts:
I am preparing an intermediate GIS course for graduate and undergraduate students who have likely not been in touch with anything but ESRI software. It is an existing course in the curriculum that I will be teaching for the first time. Currently, ArcGIS is the weapon of choice for the lab sections and practical assignments.

I want to tweak the course a little to include an introduction to Open Source GIS alternatives. For now, this part of the course will only be two to four weeks (I'm thinking a kind of extended epilogue) so I won't be able to dig too deep. I hope to branch this off into a full Open Source GIS course next year, but curriculum constraints prohibit me from doing that right away.

Here's some reading I have ben considering to support the Open Source GIS part of the course, to give you an idea of what I'm (not) looking for:

The Geospatial Desktop is a book I'd love to use for a full Open Source GIS course but is too much to cover in a few weeks. The chapter 'Survey Of Desktop Mapping Software' looks like something I could use.
The Dekstop GIS book is of similar breadth but currently out of print.
The Grass Book seems too focused on GRASS. I think GRASS will put students with an ArcGIS mindset off. Also, I want them to learn to appreciate the breadth of the OS geospatial software spectrum.
'How to go from GIS novice to Pro without spending a Dime' takes a good, practical approach to delving into OS GIS, and has good links for further reading.
A white paper from OpenGeo talking about markets for geospatial software, and how that landscape is changing.

Can you suggest other articles and / or books that would be useful to ease students coming from an ESRI / ArcGIS background into appreciating and using Open Source alternatives?


texts: What is equivalent of Esri Shapefile metadata in PostgreSQL/PostGIS database?

texts:
It is possible to enter a lot of information in an ESRI Shapefile metadata in ArcCatalog.

Is there something similar with a PostgreSQL/PostGIS database?

Is there a non-ESRI tool out there that makes it as easy as ArcCatalog to edit and view metadata of PostGIS tables?

What about ArcMap, QGIS, etc integration? It is very nice and convenient to search for a layer in ArcMap. Just search, drag and drop.


texts: Upgrading PostgreSQL with broken PostGIS installation ("$libdir/postgis-2.3" error)?

texts:
As there are two ways to solve this question, I asked for another approach here.

Whenever a major release of PostgreSQL appears, homebrew will happily install it right away, unless one pins a certain version or installs a specific version from Homebrew/core.

So far so good, but the problem is that PostGIS also keeps updating and eventually will be incompatible with pinned or Homebrew/core versions of PostgreSQL. Again, one can pin PostGIS, yet Homebrew/core does not offer previous versions. From time to time I do brew cleanup and somehow I am now left with a broken PostGIS database. Ok, so "time to jump to PostgreSQL 10.1", I though, but pg_upgrade now gets stuck with the error

  ERROR: could not access file "$libdir/postgis-2.3": No such file or directory

That means I can't back nor forth: I cannot figure out how to get PostGIS working with postgresql@9.6 and I cannot pg_upgrade to the new version of PostgreSQL (10.1).

Therefore my question can be answered in two ways. Here I am asking:

How to upgrade PostgeSQL with a broken PostGIS installation ("$libdir/postgis-2.3" error)?

The second approach is described at Installing PostGIS on older homebrew versions of PostgreSQL (e.g. postgresql@9.6)?

What I tried so far:

DROP EXTENSION postgis; does not work because the geo-spatial tables depend on it. I could DROP CASCADE but this would probably delete my tables and I would like to preserve them.


texts:
Id like to add hyperlinks to the layers in my layer control.
Leaflet site

So when you open the layer control.

Id like to be able to add a hyperlink to Basin, shoreline etc. This is my coded.

    //create basemaps variable to call within the control
        var baseMaps = {
            "Open Street Map": osm,
            "ESRI Imagery": Esri_WorldImagery,
            "OpenTopoMap": OpenTopoMap,
            "Bing Imagery": bingLayer
        };

        //create overlays variable to call within the control. Calls layer groups
        var overLays = {
            "Basin": aShp,
            "Rivers":bShp,
            "Sub-basins":cShp,
            "Shoreline":dShp,
            "Country Boundaries":eShp,
            "Partner Research Institutions":fShp,
            "Bathymetry":both
        };

        //add control to map
        L.control.layers(baseMaps, overLays).addTo(map);
        bathyImage.addTo(map);
        bathyImage2.addTo(map);


texts: Using pattern styled polygons in Tilemill

texts:
I'm just learning CartoCSS with Tilemill.

I've got a really simple requirement for styling a polygon but cannot find any hint nowhere.

What is the code for styling simple polygon fills?

For example, like can be seen here:

Black/white hatching pattern


texts:
I have .dgn files working on Bentley PowerMap, these files contain LineString geometries representing city buildings distributed in several levels (Level 1 offices, level 2 parks...), and i want to migrate them to QGIS 3.0, converted to polygons.

Can anybody help me?


texts: Query for all roads and buildings in a province of Belgium_ overpass turbo

texts:
I am looking for a query to show all roads and buildings in the province Flemisch Brabant in Belgium.

Can anyone help me?


texts: QGIS create a line graph with attributes value on a section through the map

texts:
I wonder if I can make cross sections of GIS data (other than height).
An example:

I have a shapefile of a city with building plots and attribute values of the plot size. I want to choose a random line on my map (for example: from the centre of the city to the suburb) and consider it as a cross section. Can I create a line graph that shows me the evolution of the plot size along this line?

Or a list of the shapes and their attributes ordered along the line?


texts:
When the source band having pixel values in the range of 52 and 207 is warped using the following command:

gdal_warp -t_srs wgs84 -r near in_file out_file

The resulted band has the pixel values in the range of 0 and 207. How can I make the data range unchanged during warp?


texts: Company information (address etc .) for QGIS

texts:
In order to make a donation for developing QGIS, I have to submit to my company department all "QGIS" company (company name, address, telephone, commercial registry number ...) information that will be printed on the invoice/receipt.
I search https://www.qgis.org but I can't find any information or email address to ask to.


texts: Comparing/Correlating NDVI Rasters with other index Raster

texts:
So, I've created two rasters. One for NDVI and one for another index for surface water from the same orthorectified World-View 2 image. I can visually see how the indexes are inversely proportional. Areas of higher readings in one raster correlate with areas of low readings in the other index and vis-versa.

I'm doing the processing through ArcGIS 10.2.

How would I show this to be the case in a quantifiable way?


texts: Deleting feature class after importing into feature dataset?

texts:
Here are my actions so far:

create gdb
import feature classes
create new feature dataset
import feature classes

So now I have my feature datasets with all my layers as well as the original feature classes. Is it safe to delete the duplicate feature classes?


texts: Values in attribute table not showing up

texts: I have values in an attribute table I want to label my layer with. But when I label it, a totally different set of values show up. In the attribute table, my values are between 0 and 200. But the labels have values as extreme as 300,000. What's going on? Thanks!


texts: Using Flex to display single rasters from raster catalog in ArcSDE/ArcGIS server?

texts:
One of the front end develoeprs here is having issues displaying, or working out how to display a single raster held in a raster catalog in ArcSDE using the Flex 2.2 API

the Raster Catalog holds about 30 images and he wants to search for one of these images, based on either the name of the raster, or some attribution tied to each raster, but can't figure out how.

Has anyone done this before?


texts: Selecting raster cells by points using ArcGIS for Desktop?

texts:
I'm running a spectral mixture analysis project. For my ground truth data, I've created a random points layer, and need to extract the 30x30m LandSat pixels in which the points lay. After that, I'll convert them to polygons.

So, how do I extract JUST those cells with points on them?

I'm using ArcMap 10.3.

To clarify, I'll only be using their shape, location, and orientation when converted to polygons. Their values aren't important (i.e., "Extract Values by Points" is not applicable here).

The "Extract by Points" tool DOES work, but I'd like to figure out a way to add my entire list of points, or just select the point layer, instead of manually copying and pasting all of the X and Y coordinates into the tool and adding them individually.

I'm guessing this is an ideal situation for using ModelBuilder. I have my Input Raster, a .txt file (should I use another format?). I'm stumped from here:


texts: How to transform EPSG:3857 to tile pixel coordinates at zoom factor 0?

texts:
I'm creating a Leaflet map with an WebGL overlay (using L.CanvasOverlay) which renders a street network on the GPU.

Let's assume I have a point (1497143.62, 6896952.18) in EPSG:3857 and want to transform it to a pixel coordinate on a 256x256 raster, similar to a google maps tile at zoom level 0 (compare: this and that).

I can convert the same point (52.51904, 13.44907) to EPSG:4326 and use the following JavaScript code to do the reprojection to the 256x256 tile raster:

/**
 * Converts Lat/Lon to tile pixel X/Y at zoom level 0
 * for 256x256 tile size and inverts y coordinates.
 *
 * @param {float} lat Latitude
 * @param {float} lon Longitued
 * @return Point with x and y corrdinates.
 */
function latLonToPixels(lat, lon) {
    var sinLat = Math.sin(lat * Math.PI / 180.0);
    var pixelX = ((lon + 180) / 360) * 256;
    var pixelY = (0.5 - Math.log((1 + sinLat) / (1 - sinLat)) / (Math.PI * 4)) * 256;
    return { x: pixelX, y: pixelY };
}

... which leaves me with the desired point (137.56378, 83.94214) on the raster.

So, using leaflet I could transform my point from 3857 to 4326 and afterwards calculate the pixel coordinates using the formula above. But as I understand the forumlas, this would require two spherical reprojections which is jus too much computing for large geometries.

I suspect there is a much more simple and efficient approach for this transformation. But I couldnt figure it out yet. Any idea?

In addition, searching the internet I can't figure out what the EPSG:3857 coordinates actually denote. This and that answer try to explain it a bit but it is still unclear how to programmatically solve this.


texts: Converting latitude/longitude to full range referential

texts:
I have a list of latitudes and longitudes expressed as follow:

-1 <= latitude <= 0
0 <= longitude <= 1

Those coordinates represent POIs on a fantasy map displayed using the Leaflet JavaScript library using the coordinate system EPSG:3857 (Web Mercator)

I need to map those coordinates to the full range of possible latitudes/longitudes (-90/+90 and -180/+180) in order to display them properly using the iOS MapKit framework.


texts: Snap points to the nearest location on a polygon

texts:
I have some PostGIS SQL that generates uniformly distributed point cluster centroids within polygons (loosely based on this tutorial). The number of point cluster centroids generated is based on the polygon area. The purpose of this is to create uniformly distributed sample points by area within forest stand polygons.
Here is an example of the output:

And here is the SQL (including sample data) to generate this:
-- Make up some data
CREATE TABLE polys(poly_id, geom) AS (
        VALUES  (1, 'POLYGON((1 1, 1 5, 4 5, 4 4, 2 4, 2 2, 4 2, 4 1, 1 1))'::GEOMETRY),
                (2, 'POLYGON((6 6, 6 10, 8 10, 9 7, 8 6, 6 6))'::GEOMETRY)
    );

-- Create point clusters within each polygon
CREATE TABLE pnt_clusters AS
  SELECT  polys.poly_id,
      CASE
          WHEN ST_Area(polys.geom) >9 THEN ST_ClusterKMeans(pts.geom, 8) OVER(PARTITION BY polys.poly_id)

          ELSE ST_ClusterKMeans(pts.geom, 2) OVER(PARTITION BY polys.poly_id)

      END AS cluster_id,
        pts.geom
  FROM    polys,
          LATERAL ST_Dump(ST_GeneratePoints(polys.geom, 1000, 1)) AS pts
;

-- Create centroids from point clusters, making sure all centroids are within a polygon
CREATE TABLE centroids AS
  SELECT cluster_id, ST_PointOnSurface(ST_collect(geom)) AS geom
  FROM pnt_clusters
  GROUP BY poly_id, cluster_id;

The "C" shaped polygon on the left contains centroids that are too close to the polygon boundary for my needs. Based on feedback, there seem to be two approaches to handle this situation:

Generate a negative buffer of the polygons and generate the centroids within those negative buffer areas.

Snap the centroids that are within 0.4 units of the original polygon boundary to the closest location on a negative buffer boundary

Option 1 will not work because the density of the point cluster centroids will be too biased toward the center of the polygons. For example, if I generated point cluster centroids within the green negative buffer region, you can imagine the points would be biased too far toward the center of the original forest stand (blue) polygons:

How can I snap the point cluster centroids (red points) that are within 0.4 units of their polygon boundary (blue polygons) to the closest location on the negative buffer polygon boundary (green polygon)?


texts:
I have a huge list of layers in a project that are named in a foreign language. I want to translate and add the English title to the name.

When I right click on a layer and choose Properties -> General and change the layer name, it is visualized in the layer tree but the original file stays the same.

Is there a way to change the name of the original file from QGIS?

There are so many that searching them all in my folders would take too long.


texts: I have a situation that requires the displaying of the results by precinct of an election between two of the candidates. For each precinct I would like to have two colors filling the polygon proportional to the total vote count. Also displaying the numeric count in text. Does anybody know a strategy that will achieve this result?


texts: Geometric datatype that's smaller than 32bytes in PostGis?

texts: I'm putting together a database that will contain many billions of points and multipoints and the 32byte size of the default data type (SRID=0) is starting to be a real headache for me.  Is there another SRID that stores positive integer data in a smaller size?  I'd be very happy with an unsigned 16bit integer!


texts: qgis spatial index nearestNeighbour returns arbitrary number of elements

texts:
Why does QgsSpatialIndex.nearestNeighbour(Point, 1) returns more than 1 element some of the time?

I specifically set it to return the closest one neighbour.

What's the problem here?


texts: Accuracy of raster cell size area to real world area

texts:
I read from somewhere here that a UTM raster's pixel count ratio is 1 pixel = 1 meter. I am working on a raster that says the island is 4.03 km² in area but a quick google search says the island is 4.8 km². I need to make sure that my area calculations are accurate.
How accurate is a raster's area to the real world? Does this depend on whether the data are projected on WGS 84?


texts: I need to know how to create a web map of one QGIS workspace. Actually I want to provide a link and share it to web.


texts: Can WMTS URL placeholders be used in OpenLayers 3?

texts:
In OpenLayers 2 I am using WMTS with following URL (basemap.at):

https://maps.wien.gv.at/basemap/geolandbasemap/{Style}/{TileMatrixSet}/{TileMatrix}/{TileRow}/{TileCol}.png

Is that possible in OpenLayers 3 too?


texts: How do I set map dot colors based on a JavaScript variable?

texts:
I'm trying to dynamically set my Carto map's dot colors based on a JavaScript variable.

var queryAllRows=   "SELECT+" + colNames.join(',') + "+FROM+MY_TABLE_NAME";
queryAllRows=       queryAllRows.replace(' ','+')
var queryAllRowsURL="https://ACCOUNT_NAME.carto.com/api/v2/sql?q=" + queryAllRows + "&api_key=" + apiKey;
var sums=           [];
$.getJSON(queryAllRowsURL, function(queryResult){
    // The point of this loop is to find the highest number in `rowObj`.
    // Each `rowObj` elements looks something like this...
    // {city: "New York", score1: 10, score2: 20}
    for (var i=0; i < queryResult.rows.length; i++) {
        var rowObj= queryResult.rows[i];
        sums.push(sumOfObjElements(rowObj)); // This line sums up the numbers in the object, pushes the sum into an array
    }
    var highestScore=   Math.max.apply(Math, sums);
});

Each dot on my map is a city. I want to sum get the sum of each city's score, divide it by highestScore, and set the dot color based on that result.

How do I do this?


texts: Splitting polygon with line in ArcGIS Pro

texts: I am using ArcGIS Pro and I want to cut my polygon features with a line, as the output would be also a polygon. I have no special licences, so I am looking for an alternative way. Any ideas?


texts:
I have plain data in a postgresql table which looks like this:

╔════╦═══════╦═════════════╦══════════════╦══════════╦═════════════════╗
║ id ║ index ║   easting   ║   northing   ║ utm_zone ║ utm_zone_letter ║
╠════╬═══════╬═════════════╬══════════════╬══════════╬═════════════════╣
║  1 ║   1   ║ 233410.0000 ║ 1024701.0000 ║       38 ║ N               ║
║  1 ║   2   ║ 213310.0000 ║ 1167201.0000 ║       36 ║ N               ║
║  2 ║   1   ║ 213310.0000 ║ 1167201.0000 ║       31 ║ B               ║
║  2 ║   2   ║ 213310.0000 ║ 1167201.0000 ║       31 ║ B               ║
╚════╩═══════╩═════════════╩══════════════╩══════════╩═════════════════╝

Note: Coordinates are not real, but typed by hand

As you can see these are UTM coordinates but stored in plain db type fields.
I want to convert those last 4 columns to one of Geometry type for PostGis.
To something like this:

╔════╦══════════════════════╗
║ id ║       Geometry       ║
╠════╬══════════════════════╣
║  1 ║ 02347237427342342347 ║
║  2 ║ 27584872345646325863 ║
╚════╩══════════════════════╝

I used query similar to this:

--creates polygon with wrong SRID :(
SELECT z.id, ST_Polygon(z.geom, 3395) geom FROM (
    --adds first point of line at the end to close it (requeird by ST_Polygon)
    SELECT l.id id, ST_AddPoint(g.geom, ST_StartPoint(g.geom)) geom FROM (
        -- groups points by ID and creates open line geometry from each group
        SELECT t.id id, ST_MakeLine(t.strPoints) geom FROM (
            --creates POINT(x,y) text from coordinates
            SELECT
                g.id id,
                g.index idx,
                'POINT(' || g.easting || ' ' || g.northing || ')' strPoints
            FROM data g
            ORDER BY g.index
        ) t
        GROUP BY t.id
        HAVING count(t.strPoints) > 2
    ) g

) z
where ST_IsValid(z.geom)

But this does not use utm_zone and utm_zone_letter and points and lines are created without SRID value. Creating polygon from such values and specifing some SRID gives bad results. At same time because points and lines do not have SRID set I cannot use utmzone(geometry) DB function because it requires geometry to have SRID.

I looked at ST_Transform, ST_FromText and other functions but could not find one where I can specify UTM zone number and letter to create Geometry type value.
Creating point without specifying SRID gives me bad results - points fall on penguins in Antarctica which is wrong.

Can you tell me how to convert that data to correct Geometry values with correct SRID specified?


texts: How to extract part of string after '*', '~', '-' or '#'

texts:
I have strings like this:

234/1-1
3345#5
11/2~10
6754*15

It needed only rightmost part of string after '-', '#', '~'or '*' . I figured how to extract when there is only one sign for example  '~' with right and strpos. Can I do that with all at once with field calculator?


texts:
I am sorry if this has been asked before, I have not seen a similar  question yet.

I would like to shorten a multiple attributes/fields in a shapefile in QGIS. For  example, to shorten the values like this to A- 001-F/094-A-16 to A- 001-F.

I am assuming you can do this in field calculator.


texts: Cannot install QGIS 3 on Ubuntu 18.04 (From ubuntu repo. Not ubuntugis-unstable )

texts:
When I want to install QGIS latest ver 3.3.4 to ubuntu 18.04 it rises error bellow :

Setting up qgis-providers (1:3.4.4+28bionic) ...
/usr/lib/qgis/crssync: error while loading shared libraries: libQt5Core.so.5: cannot open shared object file: No such file or directory

dpkg: error processing package qgis-providers (--configure):
 installed qgis-providers package post-installation script subprocess returned error exit status 127

dpkg: dependency problems prevent configuration of qgis:
 qgis depends on qgis-providers (= 1:3.4.4+28bionic); however:
  Package qgis-providers is not configured yet.

dpkg: error processing package qgis (--configure):
 dependency problems - leaving unconfigured
dpkg: dependency problems prevent configuration of qgis-plugin-grass:
 qgis-plugin-grass depends on qgis (= 1:3.4.4+28bionic); however:
  Package qgis is not configured yet.

dpkg: error processing package qgis-plugin-grass (--configure):
 dependency problems - leaving unconfigured
No apport report written because the error message indicates its a followup error from a previous failure.
                          No apport report written because the error message indicates its a followup error from a previous failure.
                                                    Errors were encountered while processing:

qgis-providers

qgis

qgis-plugin-grass

E: Sub-process /usr/bin/dpkg returned an error code (1)

I have removed it several times and install again but the problem still remains.

I use repository :
https://qgis.org/ubuntu bionic main


texts: Faulty Landsat 5 and 7 NDVI/ EVI archives in Google Earth Engine

texts:
For my research, I need to use NDVI or EVI datasets in Google Earth Engine. I found Landsat 5 and 7 datasets (8 day, 32 day, anual) for both. However, no matter which combination of datasets and dates I choose, I always have the issue you can see in the image below.

Looking more closely, you can see the errors along the tiles. What could be causing this? Using the Inspector tool, I confirmed the values are also effected and it is not just a visual thing...
For now I am reverting to 250m MODIS data which does not have this issue.


texts: Shapely “parallel_offset” with distance in meters from Lat/Lon coordinates

texts:
I need to find a LineString at a distance in meters from a line with coordinates in Lat/Lon. I am trying this approach based on what i found in this question https://stackoverflow.com/questions/23697374/calculate-polygon-area-in-planar-units-e-g-square-meters-in-shapely

from mpl_toolkits.basemap import Basemap
from shapely.geometry import Point, LineString
point1 = (8.83, 45.69)
point2 = (8.86, 45.71)
distance = 10000 #meters
m = Basemap(projection='sinu',lon_0=0,resolution='c')
xi_seg = LineString([m(point1[0],point1[1]),m(point2[0],point2[1])])
lcoords_s = xi_seg.parallel_offset(distance, 'left').coords[0]
lcoords_s_4326 = Point(m(lcoords_s[0],lcoords_s[1],inverse = True))

The resulting coordinates of the Point lcoords_s_4326 are (45.75149569876954, 8.745680842264079). However, i checked on QGIS and the distance between point1 and lcoords_s_4326 is lower than 10'000 meters. Any thought on how to get a distance of 10'000 meters from one point to the other? I am working with data in every part of the planet and not just in this specific region


texts: How does QGIS 3.2.3 handle or create spatial indexes within a geopackage?

texts:
My frame of reference is:

Does QGIS automatically use spatialite spatial indexes?
This explain that QGIS uses QgsSpatiaLiteFeatureIterator for Spatialite Data Provider

Spatial indexes in GeoPackage files
This explain the use of Query that is utilizing r-tree index with sub-query:

Creating Spatial Index for Geopackage using QGIS?

This last one explain that It's possible in QGIS3 to create spatial index for Geopackage. You can click "Create Spatial index" from the Layer Properties Dialog and that's it.

My doubt when reading all this, is what is the current state of the use of spatial index within a geopackage?.

Because I can use the DBmanager to create non-spatial indexes within a geopackage per column, and they are visible inside the DBmanager information window, however when using "Create Spatial index from the Layer Properties" QGIS seems to create the index -very fast by the way-, but the geopackage is not altered with any new table related to the layer.

And I can not see information in the DBmanager that the index was created.

In fact, when creating an index within the geopackage, only the option to create non-spatial indexes is possible, the option to create a spatial index appears gray/disabled.

When one creates the index from the properties of the layer, QGIS saves it inside the project but not inside the geopackage? Is this how it works?

Why is it not possible to create spatial indexes for geopackages using the DBmanager?

I'm just trying to understand the operation of the gpkg within QGIS

use QGIS 3.2.3
Data: geopackage with 2 tables
table to nonspatial: 4 million records
spatial b table: polygons approx 40k elements


texts: How to avoid clipped trash-labels in static image generated with MapBox

texts:
I choose "Generate static image" and select settings and area. All looks correct in the map. Then I enter the generated URL in a browser and look at the image; clipped trash-labels appear that wasn't there before:

For example: lower left corner, "A. D &" appears. Mid lower left above label "H mart" the text "nford" appears.
The default Mapbox Streets-style was used.

Also left an issue on github repository mapbox-studio-classic; issues 1572. Do you know of any workaround?


texts: I need to use Grid of local cooridante system at data view in order to improve georeferencing process of old map. how to make it visible, without convert it to graphics ?


texts: Trouble storing an IMap with ArcObjects 10 and IPersistStream

texts:
I'm storing a reference to an IMap in a custom object (which is stored by a custom extension) so that it can listen to map events. When I save the custom object using IPersistStream and then reload the MXD, it appears that I now have a reference to a COPY of the original IMap. I tested this by making a simple custom extension that saves the first IMap in the document... when reopening the saved MXD, the saved IMap != the first IMap in the document.

Is this how it should be working? Should I just be saving the location of the map in the document (rather than the actual IMap) and then resetting the reference during IPeristStream.Load?

I'm using C#, ArcGIS 10.0, and VS 2008.


texts: Exporting text based PDF from ArcMap using ArcPy?

texts: I'm trying to create a report style output from ArcMap or ArcGIS Pro. The first page would be a map of a state with all locations of restaurants that are in the state. Then there would be a second page that gets exported that is a list of all of the restaurants from the first page map, sourced from a field in the data table. Before I chew away at writing a python script, does anyone know if the second page is possible from ArcGIS or at least Python?


texts: Integrating ArcGIS Online Web map to my Java script file (test.html)?

texts:
Is it possible to integrate ArcGIS Online web map to my java script API?

I check the samples but in which only web map is added..

I have few layers in .html file (Like feature,dynamic etc. and I want to integrate one web map to my application.)


texts: Select parcels from SDE layer using table list and exporting to shapefiles

texts: Using a table of county fips codes and state name, I need to select from a large parcel layer stored in SDE and export to individual shapefiles based on the fips code and state name from the table.  Can this be done with a Python script or ModelBuilder?


texts:
I'm looking for some advice on how I can simplify this very detailed shapefile to ultimately be used for the web (as topojson) as a interactive heatmap. It was originally created from a raster file (geotiff - http://ecotope.org/anthromes/v2/data/). I believe it's going to take a little more massaging before using a Douglas-Peucker or Visvalingam algorithm. I believe the best way to achieve this would be to intelligently merge features based on the attribute value each contains. I'm just not sure the best way to achieve this. I have a sneaking suspicion that "dissolve" is what I need, but QGIS just freezes anytime I try that.

I'm working with command-line gdal utils and QGIS.

I have attached some screenshots below to give you a better idea of what I'm dealing with.


texts: Empty output generated when clipping features after converting KMZ to shp

texts:
I'm trying to clip two layers in ArcMap. They share spatial features and overlap each other, and they both have the same geographic and projected coordinate systems. It seems like everything should be working, but whenever I run the clip, it generates empty output.

One of these layers ("NMW boundary") was a KMZ file that I converted to a layer. Initially it was a polyline,  but after I projected it (initially, it didn't have a projection) it turned into a regular line. I am wondering if maybe that has something to do with my issue, because there isn't anything in the attribute table of the former KMZ file-- but I don't think that should make a difference.

I've attached a picture of the map I'm working on. I'm trying to clip "NG_Roads" to "projected NMW boundary", so that roads only appear within the red outline of "NMW boundary".


texts:
I am trying to work through the mapserver tutorialhttp://mapserver.org/tutorial/index.html .

I am working on a Mac (with 10.9) and I have installed mapserver from http://www.kyngchaos.com/software/mapserver.  I have apache running and I have copied the tutorial document into apache's document root ("/Library/WebServer/Documents").  All the files are world readable so apache should be able to read them. I have also verified the permissions in httpd.conf.

When I try and access the first tutorial at
http://localhost/cgi-bin/mapserv?map=/ms4w/apps/tutorial/htdocs/example1-1.map&layer=states&mode=map I get msLoadMap(): Unable to access file. (/ms4w/apps/tutorial/htdocs/example1-1.map).

Is there a separate mapserver config file that I have missed?

And where is mapserver looking for the files?


texts: How to add addresses to line segments in road network files?

texts:
I have road network data for incorporated/unincorporated places in the US. Each place is composed of road network lines with the following attributes: street name, length, place name, State.

I would like to add start/end address attributes to the lines; however, I cannot find a data sources. Does such a source exist? Or is it possible to reverse geocode?


texts: Raster tool *Extract by Attributes* creates odd output

texts:
BACKGROUND:
At 10.8.1, I have a raster tif that contains the following statistics.  Please note the maximum value:

I want to create an output raster containing those cells whose values are greater than or equal to zero.  I use the Spatial Analyst tool, Extract by Attributes as shown here:

The tool runs without error.  However, the output tif's maximum value is now different than the input tif:

This difference in maximum value is also confirmed with the Get Raster Properties tool.
QUESTION:
I expect the output raster's maximum value to be the same as the input raster.
Why are they different?


texts: Converting raster to spatial points without loss of information in R

texts:
I am currently trying to convert a rasterlayer into a spatial points object so that I can use spTransform to reproject into a new CRS. I tried to do this using the "as()" function, but the output (called Spat_points_NSW) has a different number of features (122244) compared to the number of cells (203082) in the rasterlayer object (called ndvi_median) that was used as an input.
> ndvi_median
class      : RasterLayer
dimensions : 1122, 181, 203082  (nrow, ncol, ncell)
resolution : 0.09978817, 0.008333333  (x, y)
extent     : 140.9993, 159.0609, -37.50702, -28.15702  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +ellps=GRS80 +no_defs
source     : C:/Users/U/AppData/Local/Temp/RtmpWGYBpT/raster/r_tmp_2021-05-16_201533_15536_51645.grd
names      : layer
values     : -1966, 9009.5  (min, max)
> Spat_points_NSW<-as(ndvi_median,"SpatialPoints")
> Spat_points_NSW
class       : SpatialPoints
features    : 122244
extent      : 141.0492, 153.6225, -37.50285, -28.16119  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +ellps=GRS80 +no_defs

I also note that the extent has changed during this conversion (although slightly).
So my current understanding is that for some reason, 203082-122244=80838 cells worth of data is being lost in this conversion.
How can I prevent this dataloss?


texts: Can a high GIS Stack Exchange reputation help me get hired to a GIS firm?

texts:
I have seen this question come up before, but the previous questions are all related to Stack Overflow, which is more of programming community. I am about to graduate with a bachelor's degree in GIS and I want to hear professionals input on whether or not a high reputation can impact job searching.

Should a Stack Exchange profile be incorporated into a resume? I feel like a high reputation shows enthusiasm and knowledge in the community and should be displayed for any recruiter to view, but I am just coming out of college and just getting introduced to this large community.


texts: Projection in QGIS - Scale bar and measurement

texts:
When I load a *.csv file (delimited text) as a layer in QGIS, the scale bar shows degrees (for example 0 to 10 degrees). When I uncheck On the fly projection and change the units for feet or meters, the units changed to feet or meters but still stay to 10. My scale bar is 0 to 10 meters wich is impossible because the area covered by the points is more than a thousand kilometers... But when I use the measurements tools, the measure seems to be good. Also, when I use the buffer tool, if I write 0.15 it shows a buffer of about 10 km?

I tryed many projections type and the only one that seems to work is WGS 84 but I want something more precise like NAD 83 / Quebec Lambert because the data are in the east of Canada (In the province of Quebec).


texts: For some reason I can't ever see the "Visualize" button. It flashes briefly and then switches to "Publish". This happens regardless of browser. I've updated everything and it still happens. I'm running Windows 8.1 and Chrome.


texts: How to filter point cloud by flightline in R?

texts:
The lidR package in R provides a way to mark individual points in a point cloud with a flightline ID by using lasflightline(). From the documentation:

  Retrieve each individual flightline by attributing a number to each
  point. The function depends on the GPS time to retrieve each
  individual flightline. In a continuous dataset, once points are
  ordered by GPS time, the time between two consecutive points does not
  exceed a few milliseconds. If the time between two consecutive points
  is too long it means that the second point is from a different
  flightline. The default threshold is 30 seconds. A column
  'flightlineID' is added in the slot @data

How can one filter the resulting point cloud by the flightline ID so that overlapping points from different flightlines are removed from the filtered product? Essentially, I am looking to remove the overlapping points so as to produce a homogeneous point cloud with no scanline overlaps.


texts: Vary symbol rotation to point towards fixed point using Arcade and ArcGIS Pro?

texts:
How would I write an expression in the expression builder to rotate a marker symbol towards a fixed grid reference using Arcade and ArcGIS Pro?

It does not need to be too accurate (no haversine required) as the points are within 2km and it is a projected coordinate system (ESPG 27700, OSGB)


texts: Cropping raster file with GDAL warp C++ API

texts:
I need to crop a DEM raster file with a shapefile. This shapefile has many polygons. I need to calculate several variables after the crop operation.
The algorithm is written in C++ so I want to call GDAL C++ API.
I do not need to write the warped file out to my disk and I intend to use either VRT or MEM once the test is passed.
I already obtain the coordinates of each polygon and I passed as a polygon using the cutline argument, but the operation has no effect so far.
            OGRLinearRing poExteriorRing;
            for (i = 0; i < nPt; i++)
              {
                dX = (*iIterator).vVertex.at(i).dX;
                dY = (*iIterator).vVertex.at(i).dY;
                double dDummy1 = (dX - dX_origin) / dResolution_elevation;
                long lColumn_index = long(round(dDummy1));
                double dDummy2 = (dY_origin - dY) / dResolution_elevation;
                long lRow_index = long(round(dDummy2));
                poExteriorRing.addPoint(lColumn_index, lRow_index);
              }
            poExteriorRing.closeRings();
            OGRPolygon polygon;
            polygon.addRing(&poExteriorRing);

            GDALDataset  *poVRTDS;
            poDS_elevation->GetGeoTransform(adfGeoTransform);
            poVRTDS = poDriver->CreateCopy( "./test.tiff", poDS_elevation, FALSE, NULL, NULL, NULL );

            GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();
            ...
            psWarpOptions->hCutline = &polygon;
            GDALWarpOperation oOperation;
            oOperation.Initialize( psWarpOptions );
            oOperation.ChunkAndWarpImage( 0, 0,
                            GDALGetRasterXSize( poDS_elevation ),
                            GDALGetRasterYSize( poDS_elevation ) );
            GDALClose( poVRTDS );

I assume I need to construct the polygon vertex to image index according to:
CUTLINE: This may contain the WKT geometry for a cutline. It will be converted into a geometry by GDALWarpOperation::Initialize() and assigned to the GDALWarpOptions hCutline field. The coordinates must be expressed in source pixel/line coordinates. Note: this is different from the assumptions made for the -cutline option of the gdalwarp utility!
https://gdal.org/api/gdalwarp_cpp.html


texts: Clearing Recently used coordinate reference systems using QGIS?

texts:
After using some non usual SRC for demonstration purpose, I wish to clear my  "Recently used coordinate reference systems" list for the sake of praticity.

Is it possible?

I am using QGIS 2.10.1.


texts: gdalwarp tiff reprojection from EPSG:4326 to EPSG:3857

texts:
I have a geoTiff file with the following gdalinfo output:

Driver: GTiff/GeoTIFF
Files: original.tif
Size is 48377, 15906
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (25.668509000000000,42.104685623790893)
Pixel Size = (0.000395841574135,-0.000395841574135)
Metadata:
  AREA_OR_POINT=Area
  DataType=Generic
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  25.6685090,  42.1046856) ( 25d40' 6.63"E, 42d 6'16.87"N)
Lower Left  (  25.6685090,  35.8084295) ( 25d40' 6.63"E, 35d48'30.35"N)
Upper Right (  44.8181368,  42.1046856) ( 44d49' 5.29"E, 42d 6'16.87"N)
Lower Right (  44.8181368,  35.8084295) ( 44d49' 5.29"E, 35d48'30.35"N)
Center      (  35.2433229,  38.9565576) ( 35d14'35.96"E, 38d57'23.61"N)
Band 1 Block=128x128 Type=Byte, ColorInterp=Gray
  NoData Value=15
  Image Structure Metadata:
    NBITS=4

I want to reproject it to EPSG:3857 for native (fast) rendering on Geoserver.
This is the command I use to reproject it:

gdalwarp -t_srs EPSG:3857 -co 'TILED=YES' -co 'BLOCKXSIZE=256' -co 'BLOCKYSIZE=256' -ts 48377 15906 original.tif projected.tif

Below, you can see the original, projected and both tiffs on top of each other:

Original:

Projected:

Both:

As you can see the projected raster is not aligned with the original one.

I have read that tapcan be used to align pixels, and tried this:

gdalwarp -t_srs EPSG:3857 -co 'TILED=YES' -co 'BLOCKXSIZE=256' -co 'BLOCKYSIZE=256' -tap -tr 8 8 original.tif projected.tif

Now this outputs a tiff file with a very little shift (about 2m) which I am ok with. But the processing is very slow and the resulting file is 30GB.

Giving large tr values such as -tr 40 40, fixes the slow processing problem, but the output tiles are no longer have the same value as the original ones because of the resampling I guess.

I also want to note that, tiles from the above picture are requested using EPSG:3857 srs from Geoserver. Geoserver can project native EPSG:4326 tiles to EPSG:3857 on the fly with accurate precision, but I can not produce the same result with gdal.

EDIT

These are some screenshots for different target resolutions:

-tr 20 20

-tr 30 30

-tr 8 8

-tr 7.0832868662 7.0832868662


texts: When I am trying to download my changeset in JOSM and after zooming to data in the current layer, it is showing there is no object in the data layer.


texts:
I am trying to get the coordinates of the coastline. I have made a shapefile with points selected along the coastline, the coordinate system I have the shapefile under is WGS 1984 UTM Zone 30N. I have used the Add XY Coordinates function to get the coordinates of all the points. However the coordinates don't seem to be WGS 1984, I have no idea what they are and the function has no options other than the input shapefile. Does anyone know what units the coordinates are in and how I could get them from the output that I have?

EDIT: An example of the output may help:

Point_X     Point_Y
338,682.05  5,746,608.28
338,708.81  5,746,569.18
338,750.66  5,746,540.37
338,800.49  5,746,530.92
338,839.75  5,746,543.88
338,864.39  5,746,576.43
338,891.86  5,746,620.66


texts: script to draw a square where the raster pixels contained within are a certain type and greater than X

texts:
I have a raster containing different pixel values and I am interested in having squares drawn around areas that are at least X m^2 and containing at least Y type of pixels.

The data looks like this, a coverage with the following attributes:

OID     Value     Count
0       1         123123
1       2         3342
2       4         2342342
3       5         23423423
4       7         2342323

1) I would like to identify (and draw) squares of a certain area around places where there exists at least a certain amount of a pixel value.

2) Squares may overlap.

It would be nice to have the resulting squares output to a separate layer with each square being a unique feature with the following attribute data (number / % of pixels of each value, area, dimensions...).


texts: InsertFeature on a versioned layer [vba]

texts:
I am writing an application that takes data from an Excel Spreadsheet and loads it into an ArcSDE feature class. The Feature class is editable within ArcMap, but I am having trouble inserting a feature with the code.

My code does the following

Starts an edit session.

Public Sub StartEditing(fctarget)

Dim peditor As IEditor
Dim pID    As New UID
Dim pFeatureLayer As IFeatureLayer
Dim pDataset As IDataset
Dim pMap   As IMap
Dim pmxdoc As IMxDocument
Dim LayerCount As Integer
Set pmxdoc = Application.Document
Set pMap = pmxdoc.FocusMap
pID = "esriEditor.Editor"
Set peditor = Application.FindExtensionByCLSID(pID)

If peditor.EditState = esriStateEditing Then Exit Sub

'Start editing the workspace of teh feature layer
Set pDataset = fctarget
MsgBox "Editing: " & pDataset.BrowseName
peditor.StartEditing pDataset.Workspace

peditor.StartOperation

End Sub

Load the data (if it is a new record, the following code runs)

  If pfeatLayer Is Nothing Then    'if it does not exist - it is new

  'Insert a new record and add new attribute information
  Set pFeatInsertBuf.Shape = pFeat.Shape
  Debug.Print pFeatInsertBuf.Shape.envelope.XMax
  xField = pFeatInsertBuf.Fields.FindFieldByAliasName("ORIG_X")
  yField = pFeatInsertBuf.Fields.FindFieldByAliasName("ORIG_Y")

  pFeatInsertBuf.Value(xField) = pFeat.extent.XMax
  pFeatInsertBuf.Value(yField) = pFeat.extent.YMax

    'Update the attributes
  For C = 1 To NumberOfFields
        If IsNull(pExcelApp.Cells(r, C).Value) Then
            pFeatInsertBuf.Value(C).Value = ""
        Else
            pFeatInsertBuf.Value(C) = pExcelApp.Cells(r, C).Value
        End If
 Next C
 Debug.Print pFeatInsertBuf.Fields.FieldCount
 Debug.Print pFeatInsertBuf.Value(2)
 Debug.Print pFeatInsertCur.Fields.FieldCount
 'store the record
 pFeatInsertCur.InsertFeature pFeatInsertBuf ' #######INSERT THE FEATURE######
 Print #1, "Row number " & r & ": " & "inserted"
 Debug.Print "Row number " & r & ": " & "inserted"

When it has gone through the rest of the spreadsheet, the stop editing process is called and the editor session is closed.

Public Sub stopSaveEdits(step As String)

Dim lngCityDevRefVal As Long

Dim peditor As IEditor
Dim pID    As New UID
Dim lngSave As Integer
Dim pMap   As IMap
Dim pmxdoc As IMxDocument
' get the current map document and create and editor object
Set pmxdoc = ThisDocument
Set pMap = pmxdoc.FocusMap
pID = "esriEditor.Editor"
Set peditor = Application.FindExtensionByCLSID(pID)
'step = "Save"
If peditor.EditState = esriStateEditing Then
    'lngSave = MsgBox("Do you want to save?", vbYesNo, "Save data")
    If step = "Save" Then
        peditor.StopOperation "Saving..."
        Exit Sub
    ElseIf step = "SaveStop" Then
        peditor.StopOperation "Saving..."

        If peditor.HasEdits = True Then
            peditor.StopEditing True
            'MsgBox "Changes have been saved"

            Exit Sub
        End If
    Else
        peditor.AbortOperation
        'peditor.StopOperation "Saving..."
        peditor.StopEditing False
        MsgBox "Changes have not been saved"
    End If
End If

End Sub

The problem is that the new feature never gets inserted, even though i have specified to do it. Has anyone got any ideas on what I am doing wrong?

Note - Using ArcGIS Desktop 9.3.1, ArcSDE is 9.2

Thanks!


texts: Grouping data by region in ArcGIS Desktop?

texts:
I have a map of the UK which splits the whole country into 1km polygon squares. Each polygon contains data, and I want to try and group the squares into local authority level.

I have tried overlaying a map and using a spatial join, however this has lots of double counting since the squares on the boundary are counted on both sides.

So I was wondering if there was an alternative method for doing this?


texts: QGIS 2.16 Extract Nodes command works strange

texts: On Qgis 2.16 when I am using the extract nodes command, even when I save to file with a specific name, the result is a layer named "Nodes" with the extracted points. I have to go to the folder where I saved the .shp file and load it in the program. It is not a major problem, it is somehow annoying. Do you have the same experience?


texts: Mapbox GL JS - Adding variable dependent images to popup

texts:
I've been working with popups in Mapbox GL JS, using HTML and want to be able to use variables within the Mapbox dataset, so different popups show different images, depending on attributes in the Mapbox dataset. (see below)

map.on('click', function(e) {
var features2 = map.queryRenderedFeatures(e.point, { layers: ['stadium'] });

// if the features have no info, return nothing
if (!features2.length) {
    return;
}

var feature2 = features2[0];

// Populate the popup and set its coordinates
// based on the feature found
var popup = new mapboxgl.Popup()
    .setLngLat(feature2.geometry.coordinates)
    .setHTML('
' +
        '
' +
        '' + feature2.properties['Image'] +
        '
  Stadium: ' + feature2.properties['Address'] +"
" +
        '
 Capacity: ' + feature2.properties['Capacity'] + "
" +
        '
' + feature2.properties['Fixtures'] + "
" +
        '
' + feature2.properties['Fixtures2'] + "
" +
        '
' + feature2.properties['Fixtures3'] + "
" +
        '
' + feature2.properties['Fixtures4'] + "
" +
        '
' + feature2.properties['Fixtures5'] + "
" +
        '
' + feature2.properties['Fixtures6'] + "
" +
        '
' + feature2.properties['Fixtures7'] + "
" +'

'
        )
    .addTo(map);

Basically I have Img src HTML tag, which I want to add the different image attribute depending on what is clicked. However, I've tried several ways of doing this, but can't get it to work - the pic below shows what I currently get below (without it crashing entirely). Hope this makes sense and someone is able to help...or at least let me know if it's not possible!


texts: Update OSM data already imported in qGIS

texts:
qGIS makes it really straightforward to import data from Open Street Map (there's a good tutorial here). However, despite this process is relatively straightforward, it takes some time to download the data, convert it, select the desired attributes to import, and then make queries to save a new layer.

Since I am working on a project that will take some time to develop, I am concerned about having to repeat steps or losing data if I have to make new queries from OSM to add new data. Do you know if there's any way to automate this process as much as possible.

Let's imagine that I have already downloaded OSM data, created several layers (with their corresponding attributes) and I spend some more time doing started doing some extra work (for example a heatmap). Now Imagine that I want to do the same with a more up-to-date version of OSM which will have more features. If I am not wrong I will have to repeat everything and it will be like starting from scratch. Do you know any way to just update all the layers without having to delete them and start over again? (I was wondering if there's any way to automate the process or create "smart" layers which query an specific layer or any way to update the original database)


texts: Using ArcPy to write feature class name to field in GDB?

texts:
I have a file GDB living in a network directory with 17 or 18 feature classes.  I already created a field I want to populate with the file name of the feature class.  Whenever I run my code, it returns

  NameError: name 'Washington_PETA' is not defined

This is presumably the first feature class it finds in the GDB.  File path altered because there may be sensitive data in it.

Doing this in Python 3 because we're using ArcGIS Pro.

import arcpy, os, sys

arcpy.env.overwriteOutput = True
arcpy.env.workspace = r"R:\filePathHere\feature.gdb"

petas = arcpy.ListFeatureClasses()

for fc in petas:
    print("Updating " + fc)

    inField = "PETA_Name"
    fileName = os.path.basename(fc)
    arcpy.CalculateField_management(fc, inField, fileName, "PYTHON3")

    print(fc + " completed")


texts: Adding 7 days weather forecast to QGIS project

texts: The Qweather Plugin with Yahoo Weather API gives me the real time weather info I want but I also want to add seven days weather forecast to my QGIS project. How can I implement this?


texts: Arduino compatible GPS modules with support for GAGAN SBAS and 10Hz rate support

texts:
Looking for an Arduino compatible GPS module that supports:
- GAGAN SBAS
AND
- 10Hz update rate

Any suggestions? Not able to find one that meets both requirements.


texts: Geoprocessing return not showing on map?

texts:
I have a geoprocessing service that works fine in ArcMap and ArcGIS online.

http://argis.ualr-iea.org/arcgis/rest/services/geoServices/SiteFinder/GPServer/FindSites.

I am trying to integrate it into an esri leaflet project. http://argis.ualr-iea.org/artrans/.

The service runs and returns data, but the results are not shown on the map, there are no errors and I don't recognize the Lat / Long returned?

To run this, I have an onclick handler set up, so just set up a breakpoint on line 504 (photo below) and click anywhere on the map to run the GP service.


texts: Several QGIS layers sharing single shapefile

texts: Is it possible for several QGIS layers to share a single attribute file, but each layer containing a different set of features from that one shapefile based on the value of an attribute?


texts:
I have created a SQL view from a large MultiPolygon table, and am attempting to serve it as a WFS layer via Geoserver.

While testing the layer, I've noticed that some of the features are missing. This is evident when loading the WFS layer, or also doing a SQL query on the view for the missing features (the query returns no results).

I've added 2 screenshots showing the original layer (grey polygons), and the view layer with missing features (orange polygons).

Has anyone come across a similar issue, and a fix?

UPDATE

The source table table_a has a geom field as geometry(MULTIPOLYGON, 4326). It contains about 3 million records

The view queries table_a, and does a spatial join on a couple of other tables as well, like so:

SELECT
    a.feature_id,
    a.geom,
    a.name,
    b.label,
    c.name
FROM table_a AS a
    JOIN admin_boundaries AS b ON st_intersects(a.geom, b.geom)
    JOIN places AS c ON st_intersects(a.geom, c.geom);

I'm not using maxFeatures in my query, and Geoserver is set to maxFeatures = 0

SOLVED

As suggested by @alpha-beta-soup I was using a series of INNER JOINs, and some of my features weren't intersecting, so they were missing. Switching to LEFT JOINs solved the issue.


texts: Colouring areas between (vector) contours?

texts:
Is it possible to colour areas between (vector) contours without making it a DEM first?

Joshdoe seems to have done so in Create filled contour map with GRASS but does not say how.


texts: How can I discover the Python version in QGIS?

texts:
Where can I discover the Python version installed in QGIS 2.18.1?

I would like to buy a book to introduce me at Python interface. I´ve found one that addresses the version 3.3 of Python.


texts: QGIS time manager - shapefile in edit mode?

texts:
In QGIS 2.18, shapefiles are not editable when loaded in the timeManager.
Is there a workaround?

Or do you know if this is available in more recent version?

The timestamps in the shapefile are set in seconds which by default refer to the seconds elapsed since 1/1/1970.


texts: Using File Geodatabase with Null values in Field Calculator of ArcGIS Desktop 10.0?

texts:
I want to use the field calculator to get the sum of a couple of columns. The problem is that some of them has Null Values. I am not able to find a solution.
I tried this Python functions:

stack( !basic.Blue_Thickness!, !basic.Green_Thickness!, !basic.Red_Thickness!, !basic.F1_Thickness! )

Version 1:

def stack(*args):
  myList = list(args)
  myList = [item for item in myList if (item is not None)]
  return sum(myList)

Version 2:

def stack(*args):
    return sum(filter(None, args))

Version 3:

def stack(*args):
  myList = list(args)
  for item in myList:
    if item is None:
      myList.remove(item)
  return sum(myList)

The result of all this functions is that I get always NULL back, only if there is a row with no NULL then I get a result.

Background: I have a File Geodatabase with one main table and some tables which joins this table. The calulation takes place only in the main table (source columns and write column). I am using ArcGIS 10.0.
The Source of these functions is this discussion:
Calculating field in which null values may be present?


texts: Converting Australian coordinates (EPSG:3112) for mapping in R

texts:
I'm trying to map road accidents in Australia, which comes with x and y coordinates (e.g., x = 1328459.1 and y = 1670878 for Adelaide) as well as the following information:

PROJECTION LAMBERT
UNITS METERS
DATUM GDA94 SEVEN /* GDA94 SPHEROID GRS80 PARAMETERS
-28 00 00 /* 1st standard parallel
-36 00 00 /* 2nd standard parallel
135 00 00 /* Central meridian
-32 00 00 /* Latitude of projections origin
1000000 /* False easting (meters)
2000000 /* False northin (meters)
END

Using R and this list, I tried to translate it as follows

# Projection for accident coordinates (EPSG:3112)
proj <- "+proj=lcc +lat_1=-28 +lat_2=-36 +lat_0=-32 +lon_0=135 +init=epsg:3112 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

# Code spatial points
points <- proj4::project(accidents[, c("ACCLOC_X", "ACCLOC_Y")], proj = proj, inverse = TRUE)
accidents$longitude <- points$x
accidents$latitude <- points$y
names(accidents)[names(accidents) == "longitude"] <- "LONG"
names(accidents)[names(accidents) == "latitude"] <- "LAT

before mapping, but my approach appears to be incorrect given the resulting ggplot2-map (shapefiles were obtained from Natural Earth and subsetted to Australia):

Does anyone have an idea where the error is?


texts: Selecting an attribute from a classified raster image

texts:
I am trying to select just the black attribute (classified water in the image) but I can not select anything. In the attribute table it lets me select the row, but it will not select in the map. I just want the black attribute selected (*see image).

I do not know python at all.


texts: Pulling week number from separate date column using Python and ArcPy for ArcGIS Pro?

texts:
I have a date field that contains date/time in the format DD:MM:YYYY HH:MM:SS. I want to be able to calculate the week number based on these dates in a seperate column using Python in ArcGIS Pro. I am able to get the week number based on the current date using the following code:

datetime.date.today().strftime("%W")

Unfortunately I can't translate this into pulling the week number from a separate column, rather than the current date.


texts: ArcMap 10.7 Python Label Expression - compound. expression verified but will not label. Syntax issues

texts:
I have road ID labels in a string field that need to be parsed down. I'm using the python advanced label expression. I've removed the preceding '0s' successfully as well as some known text.
I need to remove all other preceding alpha characters prior (left of) to the number but keep any alpha characters after (right of) the number. The amount of characters vary so I can not use [:3] etc.
This code snippet removes ALL alpha, not just at beginning:

This snippet shows verified but does not yield any labels:

What is your syntax suggestion?
I've also tried lstrip()[0].isalpha() can't get that to work either.


texts:
I did not find what I am looking for neither here nor in FME Website.
When I run some workflow I get that warning message (it is displayed in French I am not really sure of the accurate translation) :

  Cache functionality turn on. The workflow may run slowly because of saving features on each port and connexions

How may I turn off cache before running an heavy workflow ?

The fact is that sometimes FME stops because it has reached maximum memory.


texts: Automatically adding vertices to lines every x meters using QGIS?

texts:
I've made a vector grid using lines; the grid coverage is 100 km x 100 km and the lines are at vertical and horizontal intervals of 1000 metres. Each line in the grid is defined by two vertices, one at each end, and each line has an ordinate attribute for labelling.

My wish is to to add vertices to the lines, so that each line will have vertices at intervals of (say) 500 metres. Does anyone know how to do this?

The background to this is that the grid will be displayed on OpenLayers and given the strangeness of Google Mercator the lines need to be pegged to the ground with vertices at frequent intervals. One other obvious use for the ability to add vertices to grid lines would also be to make nicely curved lat/lon graticules for projected maps.


texts:
I'm trying to clip a raster using gdal.Translate but i get this error message:

Traceback (most recent call last):
  File "C:/Users/Mouqsit/Desktop/test.py", line 6, in <module>
    ds = gdal.Translate("C:\\wamp\\www\\Donnees_CRTS\\irrisat_dessimination-smiej\\new.tif", ds, projWin = [681483.831264, 3558183.87652, 747933.831264, 3607833.87652])
AttributeError: 'module' object has no attribute 'Translate'

This is my code:

from osgeo import gdal
import sys ,os
output_dir='C:\\wamp\\www\\Donnees_CRTS\\irrisat_dessimination-smiej'
etp = os.path.join(output_dir,"L8_L8_ETpot_24_30m_2016_223.tif")
ds = gdal.Open(etp)
ds = gdal.Translate("C:\\wamp\\www\\Donnees_CRTS\\irrisat_dessimination-
smiej\\new.tif", ds, projWin = [681483.831264, 3558183.87652, 747933.831264,
 3607833.87652])
ds = None

Have you any idea about what is happening right here?


texts: Trying to find the Qt5 version of QVBoxLayout

texts:
Migrating a plugin and started getting the error AttributeError: module 'qgis.PyQt.QtGui' has no attribute 'QVBoxLayout' I cannot find this one in the API changes I tied using qgis.PyQt5.QtGui but that didn't work.
I thin it is usage or syntax according to this
https://doc.qt.io/qt-5/qvboxlayout.html the in Qt5 the QVBoxLayout usage is a bit different. I have self.verticalLayout = QtGui.QVBoxLayout(Configuration) I am thinking Qt5 has a completely different way of setting ui's and there is a lot of ui settings in thsi; is there a converter out there for Qt4 to Qt5? Qgis2to3 didn't change any of this code.


texts: Simple projection for Cartesian plots on large-scale maps

texts:
For large-scale maps (say 1:10,000) the amount of distortion due to projection is relatively minor for my purposes.

What would be a simple way to convert latitude-longitude values to Cartesian x-y coordinates that would not require complex projection calculations to be done?

From my research so far, it would seem that the simplest option is the computation of the Mercator projection, for example as described in this answer to a Stack Overflow question.

Is that as easy as it is going to get, or is there an easier approach?


texts: Writing to Garmin IMG format from FME or other tools

texts:
At my workplace we manage protected areas. The maps on their GPSes do not update often enough. Typically enough they are the ones setting up signs for the general public.

Our rangers would like to have access to the newest protected areas on their GPS. They are predominantly using Garmin. The protected areas are represented by relativle simple polygons. Styling could be very simplified; one color per category (national park/landscape protected area/etc) and a border line.

I want to establis a workflow where protected areas are exported and made available in the garmin img-format or similar. I will probably do this using FME Server, but it seems FME can not write garmin img-files

Are there open source libraries available which can write garmin img-files which can be integrated with the FME workflow??


texts: How to automatically group values in symbology based on rules

texts: I want to display a points feature giving each point a color depending on the range of values within which it falls, because I have more than 120.000 values (one for each point) and don't want one color for each point. I know I can group values in symbology by selecting and right clicking on them, but it gets a very slow process when I have so many values, the process seems even to get stuck. Is there a way to get this grouping by assigning rules? For example I would have four ranges of values (=0, 0< value <=1, 1< value <=10, value >10). I know I could also create a new field and assign range values to each row after attribute selections, but I would avoid that.


texts: How would I create an input mask for ArcGIS Desktop?

texts: Someone asked me if we can create a simple input mask for geodatabase feature classes so that editing of attributes becomes easier. I am however not sure how this works with ArcGIS 10.x Is this possible with VBA or do I need 3rd Party tools to create masks? The last time I did something similar I could only create tools and toolbars...but no window controls...How could this be done now?


texts: How to use osgeo.Warp function in python and replace bash command?

texts:
I want to write the following python code without the bash command:

    # warp coordinates of each tile - to correct the origin of the tif files
    for file in tqdm(file_list, desc='warp'):
        bashCommand = "gdalwarp -t_srs EPSG:4647 -overwrite " + file + ' ' + file[:-4] + '.tif'
        process = subprocess.Popen(bashCommand.split(), stdout=subprocess.PIPE)
        output, error = process.communicate()

I have been trying, but I can not find a working solution. Here is my attempt to solve this:

from osgeo import gdal
input_file = 'input_file_name_string'
output_file = 'output_file_name_string'
input_file = gdal.Open(input_file)
gdal.Warp(output_file, input_file, xRes=5, yRes=5)
dst = None

My ultimate goal is to resolve the error 'Warning 6: gdalbuildvrt does not support positive NS resolution.' and correct the NS resolution with the warp function.


texts: Are there any overhead thermal images of a data center in the public domain?

texts:
I wonder how realistic it is to map the server layout in data centers based on overhead images.  A data center is a large warehouse containing thousands of computer servers, and I want to estimate the number of operating servers and their physical locations.  This seems realistic to do using a combination of visual and thermal imagery because these servers consume tons of power and generate lots of waste heat.  (The world's largest data centers consume more than 100MW each.)  To do this, I would like to find as detailed as possible overhead visual and thermal images of some data centers.

For example, several data centers publish thermal maps of their locations that look something like:

Ideally, I would like to use publicly available sources to reproduce a map like this with as much detail as possible for data centers that do not publish similar maps.

My goal is similar to projects that use thermal satelite imagery to monitor the operations of uranium enrichment facilities.  See for example the paper Monitoring Large Scale Enrichment Facilities using Thermal Imagery from Comercial Satelites from 2001.  This project used publicly available LANDSAT-5/7 data, which has thermal resolution of 120m and 60m respectively.

For reference, here is a google maps image of a google data center: https://www.google.com/maps/@33.0649146,-80.043408,16z/data=!3m1!1e3.  Is there a website online similar to google maps that can easily display a thermal version of these images? The largest building in the complex is approximately 500mx100m, so would only occupy about 4 pixels on the LANDSAT-5 images, and 16 pixels on the LANDSAT-7.    My goal would require significantly higher resolutions.  Are higher resolutions currently commercially available?  Also, this link describes airborne therm al camera systems that seem to have sufficiently high resolution for this task.  How would I find if similar airborne based images are publicly available?


texts: Measuring polyline using either QGIS or ArcGIS Desktop?

texts:
I've collected a bunch of GPS polylines (WGS84).

What is the easiest way to measure the length of these lines in KM using QGIS or ArcGIS Desktop?


texts: landsat NDVI time series plot in R, based on point or pixel

texts:
I am using R language and I would like to create a NDVI time series plot, NDVI vs time, based on irregular intervaled landsat images. I have found a method to create a raster stack and plot the average NDVI of the total area over time. Average NDVI is okay, but i would like to plot for a selected point (pixel based or xy coordinate) within the area.

Any suggestion? I am very new to R so a step by step explanation would be really helpful


texts:
I am trying to change the text in my map layout for numerous map series. I wanted to do this via Python but am running into an issue with ArcGIS Pro's new format. I used the code from other sites but for some reason I get the error
title.text = "A"
'list' object has not attribute text

This is my code
import arcpy

aprx = arcpy.mp.ArcGISProject(r'J:\238\07_MS4\GIS\NorthHaven_MS4_GIS\NorthHaven_MS4_GIS1.aprx')
layout = aprx.listLayouts('Layout1')[0]
title = layout.listElements('TEXT_ELEMENT', 'Title')
title.text = "A"
aprx.save()
del aprx


texts:
I have a raster which I want to reproject from GCS_WGS_1984 to GCS_NorthAmerican_1983. I am using the following code:

    try:
        spatial_ref = arcpy.CreateSpatialReference_management("GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]];-400 -400 1000000000;-100000 10000;-100000 10000;8.98315284119521E-09;0.001;0.001;IsHighPrecision", "", "", "", "", "", "0")
    except:
        print arcpy.GetMessages()

    try:
        arcpy.Project_management(o_dir+int_ras,o_dir+geo_int_ras,spatial_ref)
    except:
        print arcpy.GetMessages()

It gives the error:
ERROR 000732: Input Dataset or Feature Class: Dataset does not exist or not supporte

However, the input datasets most defn exist. Is there some issue with how I am reprojecting?


texts: Longitude latitude for each pixel of Tiff Image

texts:
I have a project to do and at the moment i have 2 tif images and their geografic data(longitude latitude) of their upper left pixel.What i want to do is get the common geografic place of the 2 images...so i want to get lonitude,latitude so that i can only keep the parts of the images that show the same geografic place.

The 2 images show deferent information for that place...for example the first one i a simple image and the second one shows the humidity level.
I want first to find the common place for each picture so that later i make the pictures to have the same width and length so that i can go back and forth on these 2 pictures looking at the same place but different information.

So, is my idea of finding the common place by each pixel'l geodata a good idea?

if yes...how would i do that?

if no...how should i do it?


texts:
I'm trying to run a script using python. The script requires some qgis functionality. I believe I have set the path, pythonpath, etc. correctly.
The OS is Windows 7 64-bit with qgis located at the non-standard location Z:\Program Files\QGIS 2.18.

I receive

"Z:\Program Files\QGIS 2.18\bin\python.exe" "Z:/Locating Incident
Response/Solver/pyspatialopt-master/pyspatialopt-master/examples
/binary_mclp_qgis.py"

Traceback (most recent call last):
File "Z:/Locating Incident Response/Solver/pyspatialopt-master
/pyspatialopt-master/examples/binary_mclp_qgis.py", line 11, in <module>
import qgis
File "Z:\Program Files\QGIS 2.18\apps\qgis\python\qgis\__init__.py", line
27, in <module>
from qgis.core import QgsFeature, QgsGeometry
File "Z:\Program Files\QGIS 2.18\apps\qgis\python\qgis\core\__init__.py",
line 30, in <module>
from qgis._core import *
ImportError: DLL load failed: The specified module could not be found.

Process finished with exit code 1

What am I doing wrong?

(I should note that the loading of qgis works fine in the Python console. Not surprisingly!)


texts: L.layerGroup equivalent in ESRI Leaflet for AGO feature service?

texts:
I'm fairly new to Leaflet and ESRI Leaflet, so bear with me. I'm creating four variables from the same feature service (AGO) because I want to separate them into four separate layers even though they come from the same service - based on different attribute values. I'm doing this by the where option in L.esri.featureLayer. However, when it comes time to bind the popup to each one of these variables, I want to avoid duplication and group each variable so I can just .bindPopup once because I want the popup to be the same in each layer. Is this possible with ESRI Leaflet?

See below for my code which may explain a bit better:

//create separate AUS TxDOT project layers by project status
var ausConstruct = L.esri.featureLayer({
    url: 'https://services.arcgis.com/KTcxiTD9dsQw4r7Z/arcgis/rest/services/TxDOT_Projects/FeatureServer/0',
    where: "DISTRICT_NAME = 'Austin' AND PRJ_STATUS = 'Construction Scheduled'",
    style: function (feature) {
        return {color: '#d7191c', weight: 4};
    }
});

var ausFinalize = L.esri.featureLayer({
    url: 'https://services.arcgis.com/KTcxiTD9dsQw4r7Z/arcgis/rest/services/TxDOT_Projects/FeatureServer/0',
    where: "DISTRICT_NAME = 'Austin' AND PRJ_STATUS = 'Finalizing for Construction'",
    style: function (feature) {
        return {color: '#fdae61', weight: 4};
    }
});

var ausDevelop = L.esri.featureLayer({
    url: 'https://services.arcgis.com/KTcxiTD9dsQw4r7Z/arcgis/rest/services/TxDOT_Projects/FeatureServer/0',
    where: "DISTRICT_NAME = 'Austin' AND PRJ_STATUS = 'Under Development'",
    style: function (feature) {
        return {color: '#80cdc1', weight: 4};
    }
});

var ausPlan = L.esri.featureLayer({
    url: 'https://services.arcgis.com/KTcxiTD9dsQw4r7Z/arcgis/rest/services/TxDOT_Projects/FeatureServer/0',
    where: "DISTRICT_NAME = 'Austin' AND PRJ_STATUS = 'Long Term Planning'",
    style: function (feature) {
        return {color: '#2c7bb6', weight: 4};
    }
});

//create group of project status variables for bindPopup (I think this is the Leaflet way but how in ESRI Leaflet?)
var txdotProjects = L.layerGroup([ausConstruct, ausFinalize, ausDevelop, ausPlan]);

//popup for txdot projects
txdotProjects.bindPopup(function (evt) {
    return L.Util.template('CSJ: {CONTROL_SECT_JOB}
HWY: {HIGHWAY_NUMBER}
COUNTY: {COUNTY_NAME}
LENGTH: {PROJ_LENGTH}
PROJECT CLASS: {PROJ_CLASS}
EST. COST: {EST_CONST_COST}
TYPE OF WORK: {TYPE_OF_WORK}
LET DATE: {DIST_LET_DATE}
BEGIN MILE PT: {BEG_MILE_POINT}
END MILE PT: {END_MILE_POINT}
FUND CATEGORY: {TPP_CATEGORY_P2}
WORK PROGRAM: {TPP_WORK_PROGRAM}
STATUS: {PRJ_STATUS}',
    evt.feature.properties);
});


texts: Orfeo Tool Box installation-QGIS_2.18.10-Windows 10-64Bits

texts:
I installed QGIS 2.18.10 and then osgeo4w as you suggested it. I went to the Processing Provider options and entered the paths as you described it. But I do not still have OTB in my processing toolbox.

Can you help me? I don't have any other idea so far.


texts: Flash more than one feature in ArcMap using ArcObjects and C#?

texts:
How to flash many features at same time in ArcMap using C# coding?

Currently I call this function in a loop,

private void FlashFeature(IActiveView _activeView, IFeature _feature)
{
    IActiveView activeView = _activeView;
    IFeatureIdentifyObj featIdentify = new FeatureIdentifyObject();
    featIdentify.Feature = _feature;
    IIdentifyObj identify = featIdentify as IIdentifyObj;
    identify.Flash(activeView.ScreenDisplay);
}

but unfortunately it flashes one feature at a time, and while flashing is happening it holds the process. As a result of that it flashes around 2 features in a second, instead of them all at once. I also want to flash thous features without changing my curent feature selection in map, if posable. The default ArcMap Find tool(form or whatever it is) can do it, and i'm trying to find a way to implement that in my tool.


texts: What is the difference between FeatureCollection (GEE made vs. loaded asset) in GEE?

texts:
I want to run an analysis over a FeatureCollection. It works perfectly, if I create the collection in GEE:

var p1 = ee.Geometry.Point([-122.64550158728298,38.44643840127136]);
var p2 = ee.Geometry.Point([-123.06960473413164,38.54089331487732]);
var p3 = ee.Geometry.Point([-122.85149523962673,38.33988067525676]);
var asset_points = ee.FeatureCollection(ee.List([ee.Feature(p1),ee.Feature(p2),ee.Feature(p3)]));

However, if I load a shapefile from my asset it fails:

var asset_points = ee.FeatureCollection("users/slisovski/AssetPoints");

The problem is making a Feature in L76:80:

ee.Feature(geom, {date:   image.date(),
                  name:   featname,
                  fitted: valF,
                  orig:   valO});

With the error: Error in map(ID=2_0): String: Parameter 'input' is required.


texts:
I have written a number of standalone Tools for ArcMap.

These tools are all .py files and are run outside of ArcMap. At the moment it has only been me running these tools, but I have been asked to share these with other users.

I am happy to share these, I am however not entirely happy in sharing the lessons and tricks of the trade I have incorporated into the flow of these tools.

Can I lock these files from viewing by people who shouldn’t see the text within the .py file, but still keep the script accessible to the windows task manager or a single run “double-click”?

As an aside, I would prefer to not embed these in a Toolbox, as I want to run these outside of ArcMap /ArcGIS Pro.


texts:
I have a lot of Photos with  lat, long, and altitude information.
When I start the import, I get an error message.

Do I need to prepare a  blank point layer before importing?
Is there another python plugin that is required?

I have had some Gis training, but am not a programmer.
I use win7 and Lisboa.


texts: Why can't I load more than 12 elevation files into NVIZ?

texts:
I have a large set of rasters in a GRASS mapset, which I would like to use as frames in an animation, visualizing the raster data as elevation. NVIZ is an obvious tool for this, but I'm having trouble getting started. Here's my basic approach so far. I'm starting with grey scale rasters that are actually RGB +alpha where R=G=B. I load one color channel from each raster into nviz as follows:

nviz elevation=file1.tif.blue@default,file2.tif.blue@default,...

This starts up NVIZ as expected for <12 files. With more than 12 files, the first 12 load successfully and then the process fails like this (tail end of log only):

Loading raster map <mean-L1-A1-000011_clipped.tif.blue@default>...
 100%
Translating colors from raster map
<mean-L1-A1-000011_clipped.tif.blue@default>...
 100%
WARNING: Loading data failed
WARNING: Loading data failed
WARNING: Loading data failed
zsh: segmentation fault  nviz

Do you have any idea why this is and how I can fix it? I'll need to load several hundred file at least, to pull off the animation.

EDIT: In the interest of full-disclosure, and because I might be goofing with shell commands, here's what I'm actually entering at the GRASS command line. My original post was intended as a minimal example.

nviz elevation=`cat /Users/gregory/SimData/blue_channel.txt tr '\n' ','`

Where the file blue_channel.txt contains 659 lines has the form:

mean-L1-A1-000000_clipped.tif.blue@default
mean-L1-A1-000001_clipped.tif.blue@default
mean-L1-A1-000002_clipped.tif.blue@default

In order to test cases with different numbers of files as input for the elevation parameter, I'm varying the number of files passed to head in the syntax:

nviz elevation=cat /Users/gregory/SimData/blue_channel.txt | head -10 | tr '\n' ',' | sed -e 's/,$//g'

I've also tested that the problem isn't with specific files (seemed unlikely). For example, nviz does not crash on files 1-10, or 11-20, but it does crash on files 1-20.


texts: Rededge Dual extreme pixel values from Agisoft Metashape orthomosaic

texts:
I have produced an orthomosaic in Agisoft Metashape from 2010 images using the RedEdge Dual MX 10-band multispectral sensor. The ortho was exported as a geotiff and opened in QGIS (ver 3.6.3). All 10 bands were recognized, however when I run the raster histogram calculator, all 10 bands have a spike in pixel values of 65535, with each band having a frequency count of 269.
I have attached an image of the full histogram and another zoomed in on the extreme values. Has anyone experienced this? Is there a way to remove these in Metashape prior to export? Or are there suggestions for dealing with the extreme values in QGIS?


texts: 'Generate Points Along Line' Tool in ArcGIS Pro: Measured Distance Not Matching Up

texts:
I have a line feature that, when using 'Calculate Geometry' in ArcGIS Pro, shows me it is approximately 7,000 US feet long.

I want to generate points every 10 feet on the line, and input that parameter to the 'Generate Points Along Line' tool as follows:

This means I should expect about 700 points to be generated. However, I get a count of more like 680 points in the new layer that is created.

When I go to measure the points, reason becomes clear: it shows that the points are being measured every 12.06 ft. (sometimes 12.01, 12.11, etc.), instead of every 10 ft. (ignore the midpoint here):

What am I doing wrong here? I was guessing the coordinate system/projection (I'm just using a generic WGS 84 coordinate system on a WGS 84 Web Mercator Auxiliary Sphere basemap, the default for ArcGIS Pro), but being 2 feet off every 10 feet seems drastic, and the points are on the same coordinate system as the line.

Update: I've tested this on a fresh line I drew, and this is an issue with the Calculate Geometry tool not matching up with the Measure tool, not just with the Generate Points Along Line tool.


texts: Not Able to Find ESRI.ArcGIS.Client For ArcGIS Runtime .Net SDK on WPF Samples

texts: I am trying to re-produce the demonstrated sample at this link but I am not finding the ESRI.ArcGIS.Client any where in .Net Runtime SDK references or even from "Add ArcGIS Reference" wizard from Arc Object tools but still not seeing any ESRI.ArcGIS.Client


texts: Adding SQL Trigger to update field on INSERT (multiple rows)

texts:
I have an address point Feature Class that has an addressID field. I am trying to write a trigger that will in effect auto increment this field when a new address is added to the map. I cannot alter the table and I cannot rely on users to populate this field when they edit.

So far, I have written a trigger that works great as long as only 1 edit at a time is made to the address feature class.

The working trigger does a count on the number of records already in the feature class (this includes the inserted record) and then updates the addressID field of the newly inserted record with the count.

However, if there are multiple records inserted (as in the user creates several new addresses and then saves edits) all of those records are given the same addressID.

So what I need is a modified trigger that correctly handles multiple row inserts. Problem is I don't really know how to do this and after much searching and reading many posts regarding the same on stack overflow and other forums, I still can't seem to figure it out!

The trigger that works for single row inserts is below...can anyone offer suggestions that may resolve my issue?

    CREATE TRIGGER trgNewAddressPointTestInsert ON [dbo].[AddressPointTest]
    FOR INSERT
    AS
    UPDATE [dbo].[AddressPointTest]
    SET addressID = (SELECT COUNT(*) FROM [dbo].[AddressPointTest])
    WHERE OBJECTID in (SELECT OBJECTID FROM inserted)


texts: Fetching the on-click data from the website into QGIS

texts:
I would like to fetch the on-click data into QGIS and show them as the placemarks.
Basically the issue comes from i.e this website:
https://eclipse.gsfc.nasa.gov/SEgoogle/SEgoogle2001/SE2024Apr08Tgoogle.html
whereafter the on-click you are getting the local circumstances for the phenomenon.

If, for example, I take into account, when the event start (stressed on the red), I must look for exactly the same time around the shape (which is the umbra in this case), as per below:

The issue is time-consuming since the website doesn't provide the download option for these circumstances. On top of that, there is no option do download these on-click points as a.geojson format for example.
It put me into a difficult and tedious situation, where I have to literally draw out the whole shape manually, as per below:

but this is not the end, as I cannot download the on-click points, as I mentioned above. Now if I want to get this shape externally, I have to redraw it in a different place. I used to do it in the ScribbleMaps, and save it as a .kml file afterward.

Next, I could display it i.e on Google Earth.

The same thing I can do with the QGIS if I select the Google Satellite imagery from the QuickMapServices and create the shapefile layer, where I will start drawing the shape.

How do I get these on-click circumstances into QGIS as a layer?
More details of my problem are here:
http://www.mkrgeo-blog.com/total-solar-eclipse-vizualisation-on-google-earth-part-1/
UPDATE:
I just made it a bit quicker copying the coordinates from the pop-up window. Unfortunately the process is still far too long, as I need about 180 single points to add up on my layer:

UPDATE:
The same situation applies to the points with the same magnitude, which eventually create the long-distance line, as per in the image below:

UPDATE II:
The answer below is definitely on the right track. I am just asking - is it possible to scratch the pop-up content? And make it populated in data attribute table?
How can I instantly remove the points, which I picked up wrongly?


texts: Planet Labs API: Landsat HTTP request returns duplicated scenes

texts:
Whenever i do an HTTP request for Landsat data with the following filter:

filter = {
'intersects': geojson.dumps(
              geojson.Polygon([[[-31.0,30.0],[-31.3,42.1],
                          [-6.1,42.1],[-6.1,30.0],[-31.3,30.0]]])),
'acquired.gte': datetime(2016, 01, 8, 10, 00, 00).isoformat(),
'acquired.lte': datetime(2016, 01, 8, 11, 00, 00).isoformat(),
'cloud_cover.estimated.gte': 0,
'cloud_cover.estimated.lte': 100
}

I get a response with three times the same scene. I.e. 'count' == 3, but all the three scenes are the same.

Is this a known issue?


texts: How to randomly stratify sampling within landcover classes?

texts:
I'm working on an accuracy assessment of a classified landsat image and I'd like input on the best possible way to generate random points across both the classified image and a ground truth image I created by processing vector data from aerial photo interpretation.

I've processed both images and converted them to vectors in R and would like to generate random points in specific landcovers across both images so I can compare the predicted vs. actual landcover value at each random site.

My methodology so far has been to compare both vectors (in their entirety) and create a confusion matrix on them but the images this would include the training pixels in the supervised classification so I'd rather not have them as part of the accuracy assessment. A smaller subset of assessment pixels is probably best but I'd like feedback on those who have done this in R in the past as to best practices.

What approach should I take to generate random pixels in each image so I can compare the landcover values?