texts: ArcObjects event when ArcMap table of contents selection is changed

texts:
I have been trying to figure this out for almost a full working day and I just can't get it working. All I want to do is fire an event when the user clicks different layers in the ArcMap table of contents.

I can get it working with, for example, the PageLayout using this:

IActiveViewEvents_Event pageLayoutEvent = ArcMap.Document.PageLayout as IActiveViewEvents_Event;

I have Googled and looked into a bunch of different threads, all revolving around the concept of 'IActiveViewEvents_Event'interfaces. I've been able to get the event when the selection changes within the map document itself... But I simply can't get an event to fire in the Table of Contents.

My best attempt looks like this:

IActiveViewEvents toc = ArcMap.Document.get_ContentsView(0) as IActiveViewEvents;
IActiveViewEvents_Event tocEvent = ((TOCCatalogView)toc).SelectedItem as IActiveViewEvents_Event;

However tocEvent is set to null. If I do not cast to IActiveViewEvents_Event, SelectedItem is not null...

If anyone could provide a full example of doing this I would be very grateful!


texts:
public void AddScalebar(IPageLayout pageLayout, IMap map)
    {

        if (pageLayout == null || map == null)
        {
            return;
        }

        IEnvelope envelope = new EnvelopeClass();
        envelope.PutCoords(0.2, 0.2, 1, 2); // Specify the location and size of the scalebar
        IUID uid = new UIDClass();
        uid.Value = "esriCarto.AlternatingScaleBar";

        // Create a Surround. Set the geometry of the MapSurroundFrame to give it a location
        // Activate it and add it to the PageLayout's graphics container
        IGraphicsContainer graphicsContainer = pageLayout as IGraphicsContainer; // Dynamic Cast
        IActiveView activeView = pageLayout as IActiveView; // Dynamic Cast
        IFrameElement frameElement = graphicsContainer.FindFrame(map);
        IMapFrame mapFrame = frameElement as IMapFrame; // Dynamic Cast
        IMapSurroundFrame mapSurroundFrame = mapFrame.CreateSurroundFrame(uid as ESRI.ArcGIS.esriSystem.UID, null); // Dynamic Cast
        IElement element = mapSurroundFrame as IElement; // Dynamic Cast
        element.Geometry = envelope;
        element.Activate(activeView.ScreenDisplay);
        graphicsContainer.AddElement(element, 0);
        IMapSurround mapSurround = mapSurroundFrame.MapSurround;

        IScaleBar markerScaleBar = ((IScaleBar)(mapSurround));
        markerScaleBar.LabelPosition = esriVertPosEnum.esriBelow;
        markerScaleBar.UseMapSettings();
        ;
    }

I used this code for add scale bar to my layout view.

    {
        string fileName = @"C:\Users\damitha.a\Desktop\Student\ARC2\Web_maps\exercise 2.mxd";
        IMapDocument mapDocument = new MapDocumentClass();
        mapDocument.Open(fileName, "");
        IPageLayout pageLayout = mapDocument.PageLayout;
        IActiveView activeView = (IActiveView)pageLayout;
        IMap map = activeView.FocusMap;

        AddScalebar(pageLayout, map);

    }

In this code I call "AddScalebar" method. it run correctly. But I didn't get any View. I'm new to GIS. Please tell me the way to get the Layout view.


texts: Expression not saving in QGIS Field Calculator?

texts:
When I try to enter the following expression to remove NULLS, it does not appear to be saving/doing anything.

if( "number" IS NULL, 0, "number")

Open Attribute Table
Open Field Calculator
Check Update Existing Field
Double-click Field and Values
Paste this expression on my applicable value
Press OK

My field still appears NULL for all rows.

Edit: I am working with a joined field.


texts: Why is an FMEException thrown when creating a feature using FMEObjects?

texts:
When I try to run the code below

try
{
    fmeSession = FMEObjects.createSession();
    feature = fmeSession.createFeature();
    return true;
}
catch (FMEException fe)
{
    System.out.println("Error Create Feature: " + fe.getFMEStackTrace());
    return false;
}

I get this message:

  Error Create Feature: (....\foundation\security\stfsecurity.h:578)

What might be the reason for this?


texts:
Im working with an application dealing with timezones and I hope someone here might be able to answer this timezone-related question.

Let us say that I am currently located in Norway.

The UTC-time is now: 07:00 I got the UTC-time from here: http://www.worldtimeserver.com/current_time_in_UTC.aspx.

And my time in Norway is 09:00.

So, UTC-time is 2 hours behind "Norway"-time.

I would expect that another way to express this difference would be that Norways timezone is +2 from UTC.

But it seems that im missunderstanding something here. Having a look here (world-map with timezones): https://en.wikipedia.org/wiki/List_of_UTC_time_offsets

We can see that Norway is +1.

So now lets make an example demonstrating my problem.

Im in Canada now which is UTC -6. I need to get some real-time data from Norway and I need it to display in Norweigan time.

I would start by get the UTC-time which is currently: 07:13

Next step would be to "convert" UTC (07:13) to Norweigan time by adding the correct time based on Norways timezone.

This is where I get confused..I know that Norways is 2 hours ahead...But looking at the timezone map I get the impression that its only 1 hour.

What is the logic behind this?


texts: Creating Walkshed using ArcGIS Network Analyst?

texts: I am trying to figure out the best and quickest way to show how far people could walk from the entrance of Disneyland in 2,5, and 10 minutes. I am trying to create a walkshed and do not have much experience in creating one.


texts:
I have put OpenStreetMap service in QGIS and I want to limit the area of zooming in and out. I want to let users to be focused just in one city and I don't want to permit them to pan all over the world (in the base map).

Can you provide me any link to see a Python example?


texts: Auto panning to feature from attribute table in ArcMap?

texts:
Is there a way to auto pan to the next feature in a data-set with one click.

I just used Data Reviewer and when it brings back a list of attributes that need fixing it opens up a "browse features window."

This window allows you to simply click the next arrow to pan to the next broken feature.

In ArcMap without Data Reviewer I find that I have to manually look at each entry for QA/QC i.e. right click, pan to feature, fix problem, lose place, etc.


texts:
I want to achieve the following:

Create an OpenLayers map, that runs totally in a local CRS.
The basic layer would be a raster tile layer.
The tiles will be generated by us once for each project (based on TIFFs or other images, provided by our customers) and stored in the file system (the area is very small (~ 10x10 km) and there will be none or not many redundant tiles).
Make the tile generation process as easy, as possible.

We can assume, that the images, provided by our customers have the right orientation and it should be okay, just to say: The upper left pixel of the map is x1/y1 and the lower right pixel is x2/y2 to have our local CRS. Projections are irrelevant for us here.

My question is: Is it possible in OpenLayers, to create a base layer, define the URL of the tiles and say: "the upper left corner of the first tile (0/0) is at x: 5271 and y: 12 in my local coordinate system and another point at tile XYZ is at x: 5427 and y: 634".

Or do I have to bring the tiles in a predefined format? Or do I have to create a custom projection definition or something like that?


texts:
I have a points layer with random points located along a route (polyline layer). I was wondering if there is a simple way to give a sequence number to each point in the order they show up along the route.

I am using QGIS 3.8, but I am very green and struggle with the lingo. I did install the LRS plugin hoping it might be a straight forward solution for me, but I can't make sense of it.


texts: ArcGIS Data Reviewer installation problem

texts:
I have untypical problem connected with Data Reviewer extension installation. On new computer, with fresh Windows 7 installation there is no possibility to install Data Reviewer extension in correct way:

Installation is finishing with the information "ArcGIS Data Reviewer 10 has been successfully installed." (no information about any installation errors)
Extension is not visible in ArcGIS Customize->Extensions window, Data Reviewer toolbar don't exist as ArcGIS toolbar.

Installed apps:

ArcGIS 10.0 (SP5)
Data Reviewer 10 (SP5 + Patch2)
Microsoft .NET Framework 4.6


texts:
I've created a simple model like this:

class AddPoint(models.Model):
    name = models.CharField(max_length=100)
    geom = geomodels.PointField(srid=3857)

    def __str__(self):
        return self.name

With this view I export the db contents as GeoJSON:

from django.core.serializers import serialize

def addPointGeoJson(request):
    vectors = serialize(
        'geojson',
        AddPoint.objects.all(),
        )
    content_type = 'json'
    return HttpResponse(vectors, content_type)

When I see the exported GeoJSON the SRID is not 3857 but 4326. Why happen this?

{"type": "FeatureCollection", "crs": {"type": "name", "properties": {"name": "EPSG:4326"}}, "features": [{"type": "Feature", "properties": {"name": "Napoli", "publishing_date": "2019-07-27T18:12:36.831", "updating_date": "2019-07-27T18:12:36.843", "timestamp": "2019-07-27T18:12:36.843", "pk": "30"}, "geometry": {"type": "Point", "coordinates": [14.291931129992006, 40.84290487729676]}}]}

I know that GeoJSON Serializer use 4326 as default SRID but
I expected it to change to 3857 as in the model.


texts: I am trying to convert some polyline coastline data into a polygon that I can use to clip/subset back other data sets by this polygon.  There are some small gaps in the coastline data, which FeatureToPolygon struggles to get just how I want it.  Any advice?


texts: Locating the maximum values of a raster intersection with a vector of polygons using Python

texts:
I have extracted the statistics of a raster in the polygons of a vector layer using the python zonal_stats function. I need to identify the coordinates of the found maximum values. How can I do it using Python?

This is the Python script which I used to extract zonal statistics. I need to retrieve the location of 'preci_max' for each polygon.

import json
import numpy as np
from rasterstats import zonal_stats

def max_round(x):
    m = x.max()
    if m is np.ma.masked:
        return None
    else:
        return m

vectorbac = "bacini.json"
raster = "preci.nc"

stats = zonal_stats(
    vectorbac, raster,
    copy_properties=True, geojson_out=True, prefix="preci_",
    stats=["count"], add_stats={
        "max": max_round,
   }
)

loopy=range(0,len(stats))
for i in loopy:
    print(stats[i]['properties']['NAME'],stats[i]['properties']['preci_max'])


texts: Unable to install Python 3.4 gdal bindings on CentOS

texts:
I can easily install GDAL and Python bindings on a Windows machine, but I can not reproduce the same thing on a CentOS machine. The problem is I want to use Python 3.4 (Python 2.7 is too old and outdated for me). So, I have two versions of Python on my CentOS machine and if I do

# update-alternatives --config python

I see such a picture:

+ 1 /usr/bin/python2.7
* 2 /usr/bin/python3.4

So, Python 2.7 is now active and when I do this:

# yum install gdal-python

everything works ok, except one little fact - gdal is installed for Python 2.7. If however, I switch to Python 3.4 and run the same command, I see multiple errors in the console, like

 except OSError, e:
 SyntaxError: invalid sytax

I do not know what to do with all this. What I want is just to be able to use gdal with Python 3.4.


texts: Reset Object ID After Truncating Table In SQL Server

texts:
I am working on a program to export data from a view containing a geometry type in Sql Server 2008, to an ArcSDE table on another server which also stores data using the geometry type.  The source database is not running ArcSDE.  This will be a nightly job.

Following the documentation at the link below, I am able to generate an Object ID with every new insert:

http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/in_SQL_Server/006z00000017000000/

My problem (if it is indeed a problem) is before I copy records between the two servers, I first delete all items in target table using a SQL statement:

DELETE FROM target_table_name

However, when I insert the first record after the delete, the ObjectID starts off at the next highest value from before when I deleted the table.  For example, lets assume that the highest ObjectID in the table before I deleted it was 1000.  After I delete the contents of the table, and start to load new data, the new ObjectID the dbo.i_get_ids stored procedure returns is 1001.  Each load has around 100,000 records and since this runs nightly, the ObjectID can become a pretty big value pretty quickly.  My question is, is there a way I can reset the value?  If I use the following ArcObjects code to truncate the table the ObjectID is reset to zero:

((ITable)targetTable).DeleteSearchedRows(null);

However, if there is a way to do this without ArcObjects, it would be preferred.


texts: "Add Delimited Text Layer" plugin is enabled, but does not show up in the plugin menu

texts:
I'm using QGIS 1.7.1 on Windows 7.  I want to import a delimited text layer, but that option does not appear in the "plugins" menu.  This problem persists, despite the following steps I've taken to remedy it:

Disable then Enable "Add Delimited Text Layer" plugin
Restart QGIS/Restart computer
Uninstall QGIS and GRASS using add/remove programs, then re-install

Is there any way to completely remove all plugins, then re-install and re-enable "Add Delimited Text Layer?"

Thank you


texts:
My script calls:

import pyodbc

...but returns error:

Traceback (most recent call last):
  File "C:\MyPyScripts\extract.py", line 1, in <module>
    import pyodbc
ImportError: No module named pyodbc

Failed to execute (extract).

I'm overwhelmed and confused by everything I've read on how to install this particular python module I'm trying to use in ArcGIS 10.5.

My machine is a Windows 10, just upgraded to ArcGIS 10.5.1 (from 10.3.1) and I also have ArcGIS Pro installed. I'm trying to get my work done in 10.5.1 (not really using ArcGIS Pro for much of anything yet).

My head is swirling with variations in the python between Map and Pro. Do I have to use "PIP"? (When I try running pip install pyodbc from the cmd as administrator, it says that pip is not recognized as an external command).

My questions are:

Assuming a typical ArcGIS 10.5 installation, as well as ArcGIS Pro, what python-related files/directories/runtimes/libraries would I most likely currently have installed by default? PIP?
What prerequisites is a user like me likely to need before installing a module like pyodbc?
Do I need to do something to install PIP? Or, is that only used/needed for Pro? Can I somehow install the pyodbc module without PIP?


texts:
I am trying to create polygons between contour polylines within a defined area (for which i have both polygon and polyline shps) but I seem to be thwarted at every turn. Using QGIS 2.14, I am sure it is an easy enough operation but the actual process eludes me. I have performed the task easy enough in ArcMap 10.x, but need to reproduce in QGIS/GRASS/SAGA.

As far as I know, the Polygonize tool from the Processing toolbox is the one that I'm after. Unfortunately, it spits out the same error:

'NoneType' object has no attribute 'isMultipart'


texts: Fastest way to Join tables in python automation

texts:
My issue comes form the fact that doing Joins manually is actually relatively quick (5mins), but python automation of this process is 6hrs! Here is my code:

import arcpy

GDB_MARK_TO_MARK_VECTOR_VIEW = "Database Connections\\OLE DB Connection.odc\\GDB.MARK_TO_MARK_VECTOR_VIEW"
table = "N:\\Geodetic\\GIS\\BaseLayers\\M2Mlines\\M2Mlines.mdb\\JOINtable"
feature1 = "N:\\Geodetic\\GIS\\BaseLayers\\M2Mlines\\M2Mlines.mdb\\UNJOINED"
feature = "N:\Geodetic\GIS\BaseLayers\M2Mlines\M2Mlines.shp"

arcpy.env.overwriteOutput = True

arcpy.CopyRows_management(GDB_MARK_TO_MARK_VECTOR_VIEW,table)

arcpy.XYToLine_management(GDB_MARK_TO_MARK_VECTOR_VIEW,feature1,"FROM_REF_UTM_X_M","FROM_REF_UTM_Y_M","TO_REF_UTM_X_M","TO_REF_UTM_Y_M","GEODESIC","M2M_ID", "PROJCS['NAD_1983_UTM_Zone_15N',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]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-93.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]];-5120900 -9998100 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision")

arcpy.JoinField_management(feature1,"M2M_ID",table,"M2M_ID")

arcpy.CopyFeatures_management(feature1,feature)

exit()

The XYToLine_management does not keep the table and therefore the table needs to be reattached. Is there a way to speed this up?

MAJOR BREAK THROUGH... (For me anyway) I honestly didn't think it would work but I checked all 300,000+ lines in my data table vs "the correct way" above, here is my new code that runs in 4 minutes!:

import arcpy,os

GDB_MARK_TO_MARK_VECTOR_VIEW = "Database Connections\\OLE DB Connection.odc\\GDB.MARK_TO_MARK_VECTOR_VIEW"
table = "N:\\Geodetic\\GIS\\BaseLayers\\M2Mlines\\JOINtable.dbf"
feature1 = "N:\\Geodetic\\GIS\\BaseLayers\\M2Mlines\\M2Mlines.shp"
table1 = "N:\\Geodetic\\GIS\\BaseLayers\\M2Mlines\\M2Mlines.dbf"
table2 = "N:\\Geodetic\\GIS\\BaseLayers\\M2Mlines\\JOINtable.dbf.xml"

arcpy.env.overwriteOutput = True

arcpy.CopyRows_management(GDB_MARK_TO_MARK_VECTOR_VIEW,table)
arcpy.XYToLine_management(table,feature1,"FROM_REF_U","FROM_REF_1","TO_REF_UTM","TO_REF_U_1","GEODESIC","M2M_ID", "PROJCS['NAD_1983_UTM_Zone_15N',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]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-93.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]];-5120900 -9998100 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision")
os.remove(table1)
os.remove(table2)
os.rename(table,table1)

exit()

I could still be guided to a new answer yet. I am always looking to improve, so any thoughts are welcome.

I also didn't include my breakthrough as an answer due to the fact that I don't think it answers the real question of faster Joins, it just fixes my specific problem for today.


texts: Making print composer in QGIS that uses multiple CRSs?

texts:
I want to make a map where one of the map images is in CRS 4326, and another portion of the map is in CRS 3857.

How do I do this?


texts: Selecting max value in SQL query from Extract by Attributes?

texts:
I currently have a raster that has 3 columns in the attribute table: OBJECTID, Value, and Count. As part of ModelBuilder, I am running this raster through the Extract by Attributes tool. The Extract by Attributes tool requires an SQL Where clause. I can't seem to figure out how to write something that will simply select the highest/maximum Count value.

I have tried MAX(Count), SELECT MAX(Count), "MAX(Count)", and many others. I checked the documentation page tool, but didn't find anything that helped.

The error I am getting most of the time:


texts: Set Polyconic projection to new origin for N America?

texts:
I'm helping out with an artist driven map project that requires a modified polyconic projection of the world, centered on N America. The artist wants to render globe gores, like the image attached below, but with the center on the US (around 100w 40n). In ArcMap and ArcGIS Pro I am able to modify the standard polyconic coordinate system to re center along 100w meridian, but I can't make any changes to the latitude of origin to put N America under the n point of the coordinate system.
Tools I have at hand are ArcGIS, ArcPro, QGIS and MAPublisher.


texts: Extrusion of 2d building footprints according to heights using python

texts:
I have 2d building footprints as polygonz geometries and building heights as attribute. I want to extrude them according to height values and generate building surfaces and polyhedral building records.

For this purpose in python I copy vertices of footprint and increase z values of coordinates according to heights then add these new vertices to my vertice list. So finally I achieve 8 points that represents corners of block shaped building as you can see in the picture.

How can I triangulate using these points and create triangles that represent surfaces of building using python?

Can I use a triangle library although it uses only 2d points as arguments?


texts: Using QGIS, is there a way to add a navigation panel that shows a tiny coarse-resolution (zoomed out) version of the map, to illustrate which part of the continent that the main view-window is zoomed into? Like the inset or "locator" maps that you would find accompanying any scientific figure, providing an overview of the spatial context surrounding the more detailed map. (E.g. containing a point or rectangle that visually communicates the current extent of the primary map window.) So that you can tell "where you are" without needing to toggle the zoom level of the main window? (And ideally able to be used as a tool to pan more quickly/accurately across large expanses, like the "navigator" pane in most photo image editing software.)


texts: OpenLayers 3 - display layer WFS in EPSG:4326

texts:
I have to display a layer from a WFS service (http://wms.pcn.minambiente.it/ogc?map=/ms_ogc/wfs/Bacini_idrografici.map) using OpenLayers 3.x API (..there's a reason). The WFS service is exposed in EPSG:4326.

In similar cases I have always used the following code successfully:

var sourceWFS = new ol.source.Vector({
    loader: function (extent) {
        $.ajax('http://wms.pcn.minambiente.it/ogc?map=/ms_ogc/wfs/Bacini_idrografici.map', {
            type: 'GET',
            data: {
                service: 'WFS',
                version: '1.1.0',
                request: 'GetFeature',
                typename: 'ID.ACQUEFISICHE.BACINIIDROGRAFICI.PRINCIPALI',
                srsname: 'EPSG:4326',
                bbox: extent.join(',') + ',EPSG:4326'
            }
        }).done(loadFeatures)
        .fail(function () {
        alert("error loading vector layer");
        });
    },
    strategy: ol.loadingstrategy.bbox,
});

function loadFeatures(response) {
var formatWFS = new ol.format.WFS();
  sourceWFS.addFeatures(formatWFS.readFeatures(response));
}

var vectorLayerWFS = new ol.layer.Vector({
  source: sourceWFS,
});

But now, this layer WFS isn't displayed on the map and I do not understand why?

What's wrong?


texts: MODIS product MCD43A3 stating 16 periodicity but when downloadad I have files every 8 days

texts:
I'm downloading the MODIS derived Albedo 16-Day L3 Global 500 m resolution product. I expected to have a dataset for dates every 16 days. Instead I have datasets every 8 days as can be seen in the image I uploaded.

The info on the product (link) does state:

  Phased production strategy: Produced every 8 days with 16 days of
  acquisition (i.e., production period 2001001 includes acquisition
  between Days 001 and 016, production period 2001009 includes
  acquisition between Days 009 and 024)

So every 8 days a dataset is released which contains the previous 16 days of data? so each of these datasets have a data overlap of 8 days?


texts:
The documentation (http://docs.pgrouting.org/dev/src/ksp/doc/index.html) says that the result should be 0 for the last row. Is the same 0 than -1?

I'm new on this. Thaks for the help.


texts:
This is the 32 year trend NDVI (raster data)

I would like to remove the green grids data by using R (the green pixel value is larger than 0.03)

and this one is R code that i use but i cannot remove it ...

Could anyone please suggest me what should I do ?

regards,
faris


texts:
I am trying to update an existing map in my organization without the original .mxd or feature classes.

Its going OK, but the legend is beyond anything I've seen. There's the typical ArcMap legend, and then all this data at the bottom that is essentially records from the attribute table.

Does anyone know if this is possible or if it was maybe imported from another program?


texts: Artefacts (gaps) in tiles from GeoServer

texts:
I have created a series of GeoTiff files containing global satellite raster data in EPSG 4326 format. I serve them in GeoServer (v2.15.2) and ask for the data through WMS requests from Leaflet maps in 3 different projections: EPSG 3857 (Web Mercator), EPSG 3413 (Arctic Polar Stereographic) and EPSG 3031 (Antarctic Polar Stereographic). The problem is that the Arctic map often shows artifacts in the form of rectangular areas of missing data (see image). This happens in particular when I zoom in. In the other 2 projections the data is fine. I also include an image of a single tile, showing a missing section. The missing section is always at the edge of a tile. never in the middle.

The closest hint to a solution I found here:
Artifacts in NetCDFs served to Leaflet from Geoserver over WMS when zooming

But in my case it does not seem to be related to an x=0 phenomenon. Furthermore, that post never got a final solution.

Here is an image when unticking "Enabled advanced map projections" in the Geoserver WMS settings. It shows that the tiles with artefacts touch the antimeridian with one corner and even introduces additional, triangular artefacts.


texts: How to convert AutoCAD file to Shapefile for map from ArcGIS API 3.2 for JavaScript?

texts: I have an AutoCAD file that I want to convert to Shapefile to appear in a map by API for JavaScript in .net. What do I do?


texts: Geoserver's JSON axis order and WFS versions


texts:
Can anybody tell me how to install this patch in PostGIS?

http://trac.osgeo.org/postgis/ticket/413

http://trac.osgeo.org/postgis/changeset/7534

I need to use the ST_OffsetCurve function. but it seems it's not available in the PostGIS installation(I installed the version "postgis-pg90-setup-1.5.3-2").

Can anybody tell me how to make it work?


texts: How to add an attribute value to info window as a link? (ArcGIS JavaScript API)

texts:
I am using ArcGIS API for JavaScript and I want to use a value from a feature class attribute as a hyperlink.
For example, "template.setContent("Web Address: http://www.esri.com>${Web_Address}");
This creates an info window with "Web Address" link to ESRI's website.
The web addresses for my features (about 950 features of a polygon features class) are included in the "Link" field. How do I write the code similar to the above example but the link in my info window to point the the attribute record in my "Link" field?


texts: Changing all entries in field to be entirely lower case using ArcGIS Field Calculator?

texts:
I have a field with a bunch of species names, but for a single species there a bunch of case combinations (e.g. Cat, CAT,cat). I want to make a new field which has all the entries in all lower case. I tried using the field calculator with the script = !SPECIES!.lower() but keep getting ERROR 000539.

Anyone have any idea where I am going wrong?


texts:
I installed/updated microsoft visual c++ 2013 and microsoft directX runtime 2013. After that, whenever i open qgis i get this

The program can't start because MSVCP120.dll is missing from your computer.


texts: Export to shapefile error in QGIS from POSTGIS table

texts:
I am attempting to export a postGIS layer to a shapefile in QGIS 2.6.1 Brighton. The table has 244,674 rows. The data displays perfectly well in QGIS but when I attempt to export to shapefile I get this error:

  Feature geometry not imported (OGR error: Value 1027563200 of field aland of feature 35053 >not successfully written. Possibly due to too larger number with respect to field width)

I went to the listed feature but the aland column isn't the value listed in the error. What am I missing here? Do I need to change my column settings?


texts: Unable to display values in combo-box in ArcGIS Java Script API (ver 2.8)

texts:
I am trying to display few values in combo-box in JSAPI but I am getting an error ( I am not able to view the combo-box items).

I have attached my completed code in JS-Fiddle (only removed server name)  and I searched on Google and put some alerts in that code so according to troubleshooting I found that problem is in “Store” (line no. 248 in populateList function ) and I think ESRI might updated  version 2.8 files hence getting this error.

Things that I have tried:

Inserted few alerts.
Tested it on 3 to 4 machines  (in IE, Firefox etc)
Contacted ESRI support but they told that its not their problem its dojo problem.
Want to update 3.3 versions but I think it is not support combo-box.

My application was working fine yesterday still wondering why this is happening…

I am confused now : (

Software's: AGIS DT & AGS ver.10.1 plus I am using query layer in AGIS DT 10.1

Any help will be great!!


texts:
I am a german student and study Geoinformatik.

I have the task to store Raster data in SQLite. So use the turotial from the gaia page to Create a database and store a image (TIFF) named NaturalEarth. This was no problem. In my next step I want to add this image to qgis. I try to add a rasterlayer from my sqlite database but it doesn't work. So my question is how to load the image (raster data) from my database (sqlite) to qgis?


texts: How to add another layer origin in my installed OpenLayers plugin?

texts: Can I modify my installed Openlayers plugin files to add a new layers origin? If yes, what files?


texts: Using grass:r.to.vect in PyQGIS correctly?

texts:
I would like to convert a raster to a point vector by using grass:r.to.vect from a pyqgis script.
My script gives no error as well as no output. Now I am stuck, what I have to change.

Here my code:

input = raster.tif
output = points.shp

processing.runalg('grass:r.to.vect',{"input": input, "feature": 1,"GRASS_OUTPUT_TYPE_PARAMETER": 1, "GRASS_REGION_PARAMETER":"%f,%f,%f,%f"%(xmin, xmax, ymin, ymax), "output":output})

So far I have tried:

Thinning the input layer before converting it to a raster. Manually running r.to.rast in QGIS works fine only on the thinned version.

processing.runalg('grass:r.thin',{"input": input, "GRASS_REGION_PARAMETER":"%f,%f,%f,%f"%(xmin, xmax, ymin, ymax), "output":output})

Loading the input layer as a QGIS Layer:

fileInfo = QFileInfo(input)
baseName = fileInfo.baseName()
filePath = str(os.path.abspath(input))
rlayer = QgsRasterLayer(filePath, baseName)


texts: Error in the result layer of raster calculator

texts:
I am working with QGIS 3.16.0 and I have a raster layer (A) consisting of 2 values: 1 is where there are olive tree, 0 is where there are no olive trees:

I have the second raster (B) classified the same just for a smaller area:

I need the areas of the first raster non overlapped by the second raster, so I applied raster calculator : (A) - (B), however, I did not got the difference of these two rasters, I got this result:

Which is the same as the raster I subtracted and not the difference of the two layers.
What am I missing here?


texts: Intersect analysis with a filter on attributes

texts:
I'm trying to do an intersect analysis,
# Intersect(in_features, out_feature_class, {join_attributes}, {cluster_tolerance}, {output_type})
arcpy.Intersect_analysis(layer1, output)

but I can't see how to further refine the analysis based on attributes?
layer 1 is events, and being compared against a layer of maintenance activities. When maintenance activities are finished, they are no longer relevant. So how do I want to find intersections between the two layers, where status of maintenance <> 'complete'?


texts: QGIS "Layers menu from project" plugin doesn't open files

texts:
Has anyone out there any experience with a plugin to QGIS called Layers menu from project? This plugin gives you fast access to layers that you specify in a Qgis-project file.

It won't work on my computer. It appears to be set up correctly, but won't open any shapefiles when I click on them.

Does this make any sense?


texts: Annotation doesn't update when features are inserted via FME

texts:
What I'm trying to accomplish is take two Feature Classes - one the current FC, the second with updated data.  I want to merge the changes into the original FC so its updated.  I want to update vs. completely replace the original FC because updating preserves the annotation I have placed.  I have been able to get FME to update a Feature Class in a File Geodatabase with the ChangeDetector transformer just fine, and it places a new annotation for updated/new features.  However, it doesn't follow the label expression that is set.  For example, I have a few multi line expressions set, but only the part of the expression before "vbnewline" is populating in the annotation.  Using a versioned geodatabase is not an option for this application.

The editing behaviors have been set in ArcCatalog under the annotation tab on the annotation feature to the following:

-Create annotation when new features are added.
-Update annotation when feature's shape is modified.

If I manually add a new feature to the original feature class in ArcMap, the annotation generates automatically as expected.  I can also go in and manually create new annotation for the features modified or inserted via FME and the annotation generates as expected.  Its only appears to be happening when FME does the leg work.

I am running the latest version of FME (2013 SP1), as well as ArcMap 10.0 SP2 at the ArcInfo level.


texts:
I am opening HDF5 file like this:
import gdal
hdf_ds = gdal.Open("myhdf.h5", gdal.GA_ReadOnly)
band_ds = gdal.Open(hdf_ds.GetSubDatasets()[0][0], gdal.GA_ReadOnly)

In documentation it is said that GDALRasterBand class has a method GetNoDataValue(), but i can't get a band from it. I mean no such a method on my band_ds and when I'm trying to get a band from band_ds by GetRasterBand() it returns NoneType object. Should I do something more to get a GDALRasterBand class from my band_ds?
REF: https://gdal.org/doxygen/classGDALRasterBand.html#abf5ddfd58352f8b2f98af8ab48786d9c


texts: How do I implement a JSON file with coordinates into a var with featurecollection? (TURF)

texts:
I am trying to use the nearest function in TURF, so what I have is a JSON file with information about parkades:

{
  "type" : "FeatureCollection",
  "name" : "ubcv_parking_www",
  "features" : [
    {
      "type" : "Feature",
      "geometry" : {
        "type" : "Point",
        "coordinates" : [ -123.25556236036074, 49.26489978066352 ]
      },
      "properties" : {
        "FAC_ID" : 2220,
        "FAC_DESCRIPTION" : "Old Firehall Meters",
        "FAC_ADDRESS" : null,
        "FAC_HOURPM" : "10:00 PM - 7:00 AM",
        "FAC_HOURSDAY" : "7:00 AM - 10:00 PM",
        "FAC_RATE" : "$3.50/hr Max 4 hours",
        "FAC_RATEPM" : "No Charge",
        "FAC_RATEHOL" : "$3.50/hr Max $16",
        "FAC_VERRUS" : "5604",
        "FAC_DISABLED" : null,
        "FAC_MC" : null,
        "FAC_BIKE" : null,
        "FAC_ELEVATOR" : null,
        "FAC_EV" : null,
        "FAC_UNDERPERMIT" : null,
        "FAC_GRAD" : null,
        "FAC_REZPERMIT" : null,
        "FAC_FSPERMIT" : null,
        "FAC_BOUNDARY_DESCRIPTION" : null,
        "FAC_TIPS" : null,
        "FAC_HEIGHT" : null,
        "FAC_MAP_DISPLAY" : null,
        "FAC_MODIFY_DATE" : "20131215152336",
        "PAYMENT_LINK" : " https://www.honkmobile.com/hourly/zones/5604",
        "HIGHLIGHT_FACILITY" : "0"
      }
    },

{
      "type" : "Feature",
      "geometry" : {
        "type" : "Point",
        "coordinates" : [ -123.24379874539181, 49.257371848192236 ]
      },
      "properties" : {
        "FAC_ID" : 2234,
        "FAC_DESCRIPTION" : "East Mall Road North",
        "FAC_ADDRESS" : null,
        "FAC_HOURPM" : null,
        "FAC_HOURSDAY" : null,
        "FAC_RATE" : null,
        "FAC_RATEPM" : null,
        "FAC_RATEHOL" : null,
        "FAC_VERRUS" : null,
        "FAC_DISABLED" : null,
        "FAC_MC" : null,
        "FAC_BIKE" : null,
        "FAC_ELEVATOR" : null,
        "FAC_EV" : null,
        "FAC_UNDERPERMIT" : null,
        "FAC_GRAD" : null,
        "FAC_REZPERMIT" : null,
        "FAC_FSPERMIT" : null,
        "FAC_BOUNDARY_DESCRIPTION" : null,
        "FAC_TIPS" : null,
        "FAC_HEIGHT" : null,
        "FAC_MAP_DISPLAY" : null,
        "FAC_MODIFY_DATE" : "20130912133801",
        "PAYMENT_LINK" : " https://www.honkmobile.com/hourly/zones/",
        "HIGHLIGHT_FACILITY" : "0"
      }
    },

I am trying to use the nearest function in TURF, however, on their website tells me that I would create:

var points = turf.featureCollection([
    turf.point([28.973865, 41.011122]),
    turf.point([28.948459, 41.024204]),
    turf.point([28.938674, 41.013324])
]);

So my question is how do I instead of writing out all the coordinate points of the parkades, write a code so that this var points reads all the coordinates in the JSON file.


texts:
I'am trying to make an atlas composer to get export for each line of a specific layer. Those lines are portions or roads.

In the print composer I have 2 maps, a picture frame and 3 attribute tables

I would like to only show in those attribute tables the specific line of the feature concerned by each page of the atlas. Actualy even with the "show only features intersecting atlas feature" option I got sometimes 3 lines in those tables.

Any idea on how to deal with this?


texts:
Any way to work with ARCSI (Atmospheric and Radiometric Correction of Satellite Imagery) like any other python package?

Can it be imported in a jupyter notebook?

>>>import arcsi

or is this only available through conda console?


texts: Way to add 1km x 1km UTM grid onto online map service

texts: Is there any exist services or way to add 1km x 1km UTM grid onto online map service such as Google Maps or OpenStreetMap?


texts:
I'm working with gdal and trying to run gdal2tiles.py.
The command gdal2tiles.py -z 4-18 image.tif tile works on my mac, but it does not work when I run in a linux machine.

I found the solution, already, at gdal2tiles.py gives ERROR 6 about EPSG:900913 on fresh Ubuntu 14.04 install.

But, the point is: Why MacOS recognize EPSG 900913 and Ubuntu does not, since the both code are the same?


texts: After some hunting, it seems that if I have an a geojson file, I NEED the data and map to be on a server in order to view this in OpenLayers because of the protocol class. I saw an example where I can basically read the file as a javascript object by inserting it into my code. But there must be a way to simple throw my .json file without a server - right? I was able to do it with leaflet... I'm a newbie at this so any help would be much appreciated!


texts: Serving multiple styles openstreetmap-carto(4.x versions) and osm-bright from tile server?

texts:
I have followed this link to implement serving multiple tile styles from tile server. It says that I just need to update .xml style file path in renderd.conf file. But the openstreetmap carto recent 4.x releases requires specifying --style and --tag-transform-script(lua) during osm2pgsql data import. And if I do that, osm-bright tile serving is not working.

I tried -

sudo osm2pgsql --drop --hstore-all --slim -r pbf -C 32000 --style  /data/
/styles/osm/openstreetmap-carto-master/openstreetmap-carto.style --tag-transform-script  /data/
/styles/osm/openstreetmap-carto-master/openstreetmap-carto.lua  --flat-nodes
node.cache --number-processes 8 -d gis -U root /data/nydata/new-york-latest.osm.pbf

After this, if I run-

renderd -f -c /usr/local/etc/renderd.conf

It is showing errors for osm-bright style because osm2pgsql query contains  openstreetmap-carto.style instead of default.style(default.style is required for osm-bright).

So my question is what should be the correct osm2pgsql query for .pbf data import so that I can use same database for osm-bright and openstreetmap-carto styles both? Or What should be the steps to make a tile server that serves both openstreetmap-carto and osm-bright tiles?


texts: Image viewer or browser for large images

texts:
I have some big images. I'm looking for an easy way to browse and view them. Currently I'm using Preview.app on OS X, but it wasn't designed for these types of images.

My images are, for example:

.JP2, 42 MB, 50,000 x 30,000 pixels
a 3.4 GB TIF
something in NITF format

etc.


texts: Does a common algorithm exist to find the midpoint of a linestring by distance within a coordinate system?

texts:
Recentl I needed to find the mid-point of a linestring. I work with java/Scala so I naturally reached for the JTS library but couldn't find anything in there for doing this. I implemented an algorithm myself and am wondering whether there is an equivalent within common GIS programming packages, or if it is a specific solution to this problem that a more general algorithm could also solve. Here's a diagram of how it works.

The basic idea is to precalculate the length and divide by 2 so you know the distance to travel to the middle. Then, you traverse the linestring from one end until you have traversed more than the midpoint distance. Finally, take the vector which points from the too-far point back to the previous point and travel the distance equivalent to the difference between how far you've traveled and the total length to the center.
Here's the Scala code using JTS.
  val distanceToMiddle = gf.createLineString(coordinates.toArray).getLength / 2
  val result = findCoordinateByTraversing(coordinates)(distanceToMiddle)

  def findCoordinateByTraversing(points: List[Coordinate])(targetDistance: Double): Option[Coordinate] = {
    findPointByTraversingHelper(points, 0)(targetDistance)
  }

  @scala.annotation.tailrec
  private def findPointByTraversingHelper(points: List[Coordinate], distanceAccum: Double)(
      targetDistance: Double
  ): Option[Coordinate] = {
    points match {
      case pointA :: pointB :: points =>
        // If the total distance travelled puts us beyond the middle, find the midpoint
        val distanceTravelled = distanceAccum + new LineSegment(pointA, pointB).getLength
        distanceTravelled match {
          case d if d > targetDistance =>
            // Find the midpoint by backtracking from pointB to pointA in the direction of pointA and magnitude of the remaining distance
            val directionAngleRadians = Angle.angle(pointB, pointA)
            val remainingDistanceToMid = distanceTravelled - targetDistance
            val newX = pointB.x + remainingDistanceToMid * Math.cos(directionAngleRadians)
            val newY = pointB.y + remainingDistanceToMid * Math.sin(directionAngleRadians)
            Some(new Coordinate(newX, newY))
          // Else keep looking for the midpoint
          case _ => findPointByTraversingHelper(pointB :: points, distanceTravelled)(targetDistance)
        }
      // Single point or no points results in nothing being returned
      case _ => None
    }
  }


texts: Change polygon's edges based on another feature

texts:
I'm trying to map out fishing areas. I had a description of their western and eastern limits, and min/max operating distance from the shore, which I turned into a list of coordinates for their NW, NE, SW and SE corners. I then generated polygons shaped like parallelograms, which worked for the most part, because the coastline in my study area is usually in a straight line. However, there are some areas where it's not the case and my polygons overlap with the coast, like so:

I also have other polygons that don't overlap with the coastline, but don't parallel it and therefore don't represent the correct min/max operating distance from shore throughout their whole area.
I have a line feature of the coastline itself, and several buffer features that represent different distances from the coastline. Is there a solution where I can use these features to easily/automatically re-shape the polygons so that they are at the correct distance from shore throughout their whole area, not just at the corner vertices? This would have to identically modify both the top and bottom edge of the polygons so that they remain parallelograms.


texts: How to configure Geoserver Geowebcache to use custom projection?

texts:
I'm running Geoserver 2.3.2 with integrated GWC on my Linux computer on Tomcat. I created a custom gridset via Geoserver GUI, which I want to use for seeding, however I can't see it in grid sets combobox on gwc page.

I've read something about geowebcache.xml, but I can't find it in the directory structure neither. There is also no gwc dir in data-dir. Do I have to create these manually?

Also, when I try to seed the layer with the default epsg:4326 gridset, the task seems not to start running. The page is reloaded but no tiles are generated.

Any help appreciated.


texts: I imported a few geotagged images and created a points shape showing their locations (or more accurately, the location of the camera) with the "Geotagged Photos to Points" tools. The attribute table contains a column called direction, which contains the direction the camera was facing. Is it possible to show that direction with arrow icons?


texts: Error merging 2 point layers with same CRS

texts:
I've edited 2 point layers with the Refactor Processing tool in QGIS 2.18.5 to re-order and edit the attribute fields so both layers will merge without creating additional fields in the merge. With both point layers saved with the same CRS, the merge tool gives me the error that 'All layers must have the same geometry type!'

I've also checked the field types match up (integer, string etc) which they do.
I've tried with the raw data, before changing the order of fields in the attribute table but still no luck. I currently have no issue merging other shapefiles.

An additional observation: I've noticed when my data goes out of the field of view within the map canvas, the icon in the layers panel changes from a single dot to 3 smaller, unfilled dots (the same icon when categorising the styling of point data).


texts: SAGA-GIS:Viewshed creation from SHP instead of interactively?

texts:
I am using Saga-GIS to calculate a viewshed from a single point. I can do this with the interactive module, (terrain analysis > Lighting > visibility (single point) interactive)  but how do I do it from a shapefile?

thanks for any tips,

Robert


texts:
In QGIS I have used the OpenLayers plugin to locate a wind farm. I have then used the add feature function to place a point at the base of 7 turbines. I have saved the shapefile as an ESRI shapefile.

I am trying to open this shapefile in OpenWind but am getting back the following errors:

"failed to open file"
"Problem reading SHX file - aborting Shape.Read()"
"error loading vector layer"


texts: Finding overlapping polygons in two shapefiles and deleting them in R?

texts:
I have two shapefiles with more than 1000 polygons each, and they represent two types of agricultural fields in the same area, which means they should be different types.

Some of the polygons are partly or completely overlapped - what I want to do is to delete those polygons - from both shapefiles - that are overlapped. Are there any methods for doing this in R using pachage sf or rgdal?


texts: What is the highest obtainable resolution of FREE elevation data in the UK?

texts:
I am seeking closure with this question guys.

I am aware of various sources of free elevation data ~30+ m, and I have been searching for higher resolution FREE data for a while now. This already exists somewhat (Where to acquire elevation data for UK?), but may have changed since a year.

I'm afraid I may know the answer already, but to put an end to my misery and false hope:

What is the highest resolution, freely available (outside of academic institutions), elevation data covering the ENTIRE UK? - is 10 m available?


texts:
I want to automatically generate a series of prints following a street for planning communcation lines next to it.

Requirements:

fixed layout
fixed map scale
serial print following selected objects of a layer or a set polygon
text-chart linked with attributes of layers.
pages should overlap ~5%
prefer pdf as print output

I noticed I can manually create and print a document with multiple pages and maps using the atlas function.
I'd be fine with a screenshot showing the settings and how to set find these options.

Following image shows in a Index page how it should set up the pages.
The red lines will be drawn as an layer
Is it possible to also get the index page printed in QGis

A single Page should look like this

This is my current Atlas setting window but I don't see an option to let QGis generate an Atlas by itself with the coverage layer.


texts: Exporting landcover clasiffication from ENVI to SAGA

texts: I have basic knowledge in remote sensing and GIS. I have generated landcover map in ENVI 4.7. How can I transfer or Export/Import this into SAGA so that I perform other functions in SAGA. I want to generate the factor (R*K*LS*C*P)in RUSLE (soil loss modelling) using SAGA.


texts:
I'm trying to get an Imagemosaic of Sentinel2 data with one coverage for each band to work.

As described in the docs, I want to make use of the CoverageNameCollector property. But the docs are very vague about it and I cannot get it to work. I always end up with just one Coverage.

Beisdes the official docs, I did not find any other source for documentation or tutorial.

My indexer.properties file looks like this:

GranuleAcceptors=org.geotools.gce.imagemosaic.acceptors.HeterogeneousCRSAcceptorFactory
GranuleHandler=org.geotools.gce.imagemosaic.granulehandler.ReprojectingGranuleHandlerFactory
HeterogeneousCRS=true
MosaicCRS=EPSG\:4326
TimeAttribute=time
PropertyCollectors=TimestampFileNameExtractorSPI[timeregex](time),CRSExtractorSPI(crs)
CoverageNameCollector=org.geotools.gce.imagemosaic.namecollector.FileNameRegexNameCollectorSPI:regex=(B[0-1][0-9]).jp2$
Schema=*the_geom:Polygon,location:String,time:java.util.Date,crs:String
CanBeEmpty=true

Example Filename: T60VXJ_20170421T232711_B03.jp2

I use current development snapshot of GeoServer (2.12).


texts: How to move a layer in QGIS Toc above an existing group?

texts:
If in QGIS TOC a group already exists and I add a new layer via

QgsVectorLayer* lyr= mQGisIface->addVectorLayer(uri->uri(),Title, myProviderName);

the new layer will be inserted within the group. How can I move it outside/above the group?


texts: Creating a cumulative factor based on multiple buffers

texts:
I'm using QGIS 2.12.0 and am trying to formulate a hierarchy for the risk factor associated with footpath polygons based on their proximity to pedestrian generators and the importance factor of that generator. For example, a bus stop will have a buffer radius of 100m and a factor of 10, while a small park will have a radius of 400m and a factor of 5. When the buffers are created, it's obviously going to create several hundred overlapping buffers and I wanted each intersecting segment to inherit a cumulative factor based on the overlapping factors (I've tried to illustrate what I mean in the pic).
Heat maps using the factor work ok, but they don't take the radius into account.

Does anyone have any ideas how I might go about achieving this?


texts: Problem with a function to create a unique image starting to image collection

texts:
Starting from Global Hansen dataset and in particular lossyear map, with the following function I obtained an image collection. Each different image is given by the overlap of image for the current year and the image of the following year. These two maps are obtained following two different condition as report.

var mylist = ee.List.sequence(1, 16, 1);

var combine_img = ee.ImageCollection.fromImages(mylist.map(function(y) {
  y = ee.Number(y).toInt();
  var a = gfc.select(['lossyear']).eq(y).rename('a'); //Current year
  var aP = gfc.select(['lossyear']).eq(y.add(1)).rename('aP'); //Next year
  var b = mapB.select(y.subtract(1)).rename('b');
  var bP = mapB.select(y).rename('bP');
  var c = mapC.select(y.subtract(1)).rename('c');
  // mapC is the discriminating factor.
  // In the map of following year I would like to have pixels that result by filtering for this map.

  var b_gt_bP = a.updateMask(b.gt(bP)).updateMask(c);
  // First condition
  // Values of map b  are greater than map bP
  // .gt mean "greater than".
  // Event in the current year

  var bP_lte_b = aP.updateMask(b.lte(bP)).updateMask(c);
  // Second condition
  // Values of bP map are greater than map b
  // .lte means "less than equals"
  // Event in the following year.

  return b_gt_bP.blend(bP_lte_b).set('imageNumber', y);
  // the first map is overlap to the other map.
  // Obtain a map with pixel value equal to imageNumber
}));

How can I modify this function to obtain not an image collection but a unique image?
In this map the pixels values are equals to imageNumber?

EDIT: In my mind the pixel values have only integer number starting from 1 to 16.

I attached a picture above. The image shows the proprieties of the first image that makes up my collection. The pixels in this image have value of 1. In my final map I want that these pixels have a value of 1 because it is referred to 2001. In the same map I would like also the pixels of other image of my collection that have a values between 2 to 16 according the year.

Here the code https://code.earthengine.google.com/b5594a0817fdd8c7483ca4906abf9eb9


texts: python qgis.core import error caused by spatialite libraries

texts:
I have an error when trying to import qgis.core to python:

>>> import qgis.core
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ImportError: /usr/lib64/libspatialite.so.2: undefined symbol: GEOSBufferParams_destroy

Also DB Manager in QGIS can't work with spatialite due to this error: "/usr/local/lib64/pyton2.7/site-packages/pyspatialite/_spatialite.so: undefined symbol: GEOSSimplify"

pyspatialite seems to work fine:

>>> from pysqlite2 import test
>>> test.test()
.....................................................................................................................................................................................................................
----------------------------------------------------------------------
Ran 213 tests in 2.892s

OK

The issue arose after updating QGIS from 1.7 to 1.8, where pyspatialite (and some other packages) is no longer provided within QGIS.

I searched the web but wasn't able to find sollution.  Any idea on how to fix it?

OS: openSUSE 12.1; python 2.7.2-7.17.1; QGIS 1.8.0-7.2; spatialite 3.0.1-9.1; GEOS 3.3.5-1.1; pyspatialite 3.0.1.


texts:
There is a default chloropeth map in GeoServer:

Is there a way to style a Choropleth map based on Four_color_theorem so the adjacent sides would always be different colour?


texts:
In GeoNetwork I try to harvest a directory containing some geographic files. The harvesting seems to work fine but I can't see records on the catalog.
Is there any documentation I could follow? Or do I have to load data into a GeoServer first?


texts: Selecting lossyear polygon of Hansen data

texts:
I'm working with Hansen database in Google Earth Engine.

I would like to select polygon of two different years for example 2010 and 2011 for a certain area using lossyear map.

This is my code:

var hansen =ee.Image('UMD/hansen/global_forest_change_2017_v1_5')
.clipToCollection(area)
.select(['lossyear']);
var 10 = loss_year.eq(10);
var 11 = loss_year.eq(11);

To create a unique map I think to use the command .and() so

var map = 10.and(11);

When I display the result, the map is not visible why?


texts: How to change the aspect ratio of a map without the scale changing?

texts:
The aspect ratio of the map frame doesn't match the contents, leaving a lot of white space inside the frame, even after I've scaled the contents to fill one dimension. When I drag one edge of the frame to change the aspect ratio, the scale of the content changes. I know I can zoom with the mouse wheel, but the jumps are too big. Is there a way to set the sensitivity of the zoom? Alternatively is there a better way to change the map's aspect ratio so that the scale doesn't change?
qgis 2.01 64 bit
windows 7


texts: Adding a start and end marker to a GeoJSON linestring

texts:
I have a GeoJSON route that was created using ORSM which is as follows:

var noelev = {"type":"Feature","properties":{},"geometry":{"type":"LineString","coordinates":[[-84.477478,37.988621],[-84.477786,37.987948],etc,[-86.448045,36.984769]]}}

I add the data to the script using:

var noel = noelev;

I add this data to the map and it shows the line string

 L.geoJSON(noel, {

       }).addTo(map);

And I am trying to add a marker point to the start and end of this line

I have tried using:

   var coords = noel.features[0].geometry.coordinates[0];
   var marker = L.marker(new L.LatLng(coords[0], coords[1]), {

   }).addTo(map);

I get a TypeError: noel.feature is undefined
But this does not work. I have tried everything!


texts: Not able to load a DEM raster into Spatialite

texts:
I am using Spatialite to store a DEM raster with terrain values (GMTED2010 from Nasa). To be more precise, I am using the version: spatialite_gui-1.7.1-win-x86.zip plus rasterlite_load (rasterlite_tool-1.1g-win-x86.zip) and rasterlite_pyramid (rasterlite_pyramid-1.1g-win-x86.zip). The problem is that the DEM raster gets into Spatialite as an RGB raster, so it does not longer have the associated altitude values but 3 values (RGB) per pixel instead, which makes the appearance of a greyscale DEM but without heights!. I am using the following syntax:

rasterlite_load -d GMTED2010.sqlite -T raster -D . -i tiff

So, I thought Spatialite handled the type of raster, but apparently it only uses RGB types, doesn´t it?

What I am missing?


texts: Is it possible to use comments in QGIS Query Builder with ESRI shapefiles?

texts:
I am trying different queries in the Query Builder.
Is there a way to comment the query?
I tried with

/* */

#

--

without success :

  An error occurred when executing the query.
  The data provider said:
  OGR[3] error 1: SQL Expression Parsing Error: syntax error

EDIT :

I have made a feature request of this question: http://hub.qgis.org/issues/10914
@Bernd V. observed that expressions are evaluated differently according to the format of the vector (see the [help] button in the Query builder for more information)


texts: QGIS Draw dynamically marker's heading from field in a postgres database

texts:
General description:
I am using QGIS as quick visualization tool, where I make a "PostGIS" database layer and I see user generated events.
Table structure as follows:

Id   |   time_of_event  | position         | heading  |
int  |   timestamp      | geometry(Point ) | smallint |

Specific problem:
When I draw the position I draw it as point, I would like to draw it like so, if heading is = 0 or > 360, otherwise I would like to draw it as arrow.

I have read How to draw line pointing in compass bearing direction from point shapefile and Scaling vector field marker arrow based on attribute? but both of them are static.
And I couldn't figure out how to draw them based on a conditional, is this even possible?

P.S. I am using QGIS 3.2

Edit: so I made it display the proper direction:

However how can I make a conditional, when heading is null, display red dot?


texts: Assign values to a subset of cells of a raster

texts:
I have a data frame (3xm) with longitude, latitude and a parameter value (population density) as columns
 These points, long and lat,  correspond to the centroids of the cells of a gridded raster. The data frame has not all of the centroids but  only the ones whose cell has  the parameter value different from zero.

I have the full gridded raster in addition to the data frame

How can I create using R a new raster with the same cells of the raster that I have  where the values in the data frame are assigned to  the cells whose coordinates are in the data frame  and NaN to the other cells?

Note that the dimensions of the raster is very large


texts: ArcMap Raster Calculator SetNull not changing value?

texts:
I am having issues with the "SetNull" function within the Raster Calculator in ArcMap.  I have a landsat 8 image that I separated the bands into their own raster (only working with 6 of them), bands 2-7.  I needed to separate the bands beucase each band has a different "NoData" value. I.e. Band 1 "NoData" = 10, Band 2 "NoData" = 5, so on and so forth.  I was able to perform the following on all the bands except for band 2:

Raster Calculator
SetNull("Band1" == 10, "Band1")

Band 2 would run, but the "NoData" pixels would stay at their same value (would not change).  The only value that is changing is the stretch value.  See images below.  The red dot was painted on to show which pixel I am picking.

I've tried a couple different operations with band 2 such as trying to multiple by 0.5.  the pixel value will change then, but setting it to Null does not work.

Any suggestions?


texts: Don't add output result to map in ArcGIS

texts: Whenever I execute some GIS operation the output is automatically added to the map. Is there any option in ArcMap to turn off this setting?


texts: Is it possible to link symbol fill colors between different geometry types in QGIS

texts:
I would like to link symbol colors between two layers of different geometry types (point and polygon). I have tried to save/load style and to copy/paste but I don't think its designed to work that way across geometries.

Specifically, I have a layer of irrigated parcels with a field for the diversion structure code (ditch headgate) that waters that parcel. I also have a ditch point file with an attribute for each ditch number. More than one parcel can be irrigated by a single ditch ID.

The ditches occur at different locations across the landscape are geographically separated from the parcels.

I would like to color the ditch point fill and the label text the same as the parcel they are irrigating so it is possible to visualize where the water source for each parcel is located.

Possible in qgis to dynamically link the point layer fill colors to be identical to the random colors in the polygon layer?


texts: Isometric map tiles with Openlayers or similar

texts: I wish to show a map layer in a browser based tool such as Openlayers or Leaflet however my tiles are projected isometrically.  Has anybody done this before?


texts: Style WM(T)S layers from a QGIS server into a Cesium imagery component

texts:
I have a QGIS server docker image which is serving some aerial images (raster data) on one virtual machine.
I have an other virtual machine which holds a Cesium cartographic web application in Vue.js.
I am wondering if it's possible for example, to dynamically apply different styles to, e.g., a Cesium WebMapServiceImageryProvider which URL parameter points to the QGIS server available layers, as if I was styling these layers directly in QGIS?


texts: Find out if overlaying polyline matches polygon centerline in qgis

texts:
I have a polygon and a polyline and I need to find a scientific way of how close that polyline is to the actual polygon centerline. Is there any way of doing this in qgis? Any extensions, models, algorithms would be helpful!

Thoughts about this in my mind so far:

Calculate centerline of polygon and then calculate distance of each point to nearest point of the other polyline. Problem: Calculating Voronoi polygons in qgis is crashing on my machine.
Buffer that polyline, substract from that area the given polygon to calculate the difference area. The smaller, the better. Problem: I have no attributes in that polygon to distinguish between classes and so thin roads will give a different result than broad roads.
Do a negative buffer so that the area is actually really small and clip it with the polyline. Same problem as above.


texts: How to Get the Three Nearest Points to an Object (different layer)?

texts:
I have two point shapefiles and I want for each point in the 1st shapefile to get the 3 closest (based on air distance) points from the 2nd shapefile.

If possible, it would be great to have each of the three as a new column (as if a variable name) and the value is the distance.

Thanks a lot!

M


texts:
I have a similar problem to the user in this question from a couple of years ago:

How to Extract by Attribute Using ModelBuilder with User Input?

Except when I go to build the model, the "Values" variable does not exist in the dropdown. Is there a newer way to do this now?

I have, essentially, a parcel layer that I want the user to input the digits of each parcel they want to extract and then export that to a new layer.

In the past, I just built the SQL in Excel and copied that into the Query Builder and then exported the data from there. However, I want to create this tool that a few of my less computer-oriented coworkers can use so it eliminates third party (ie Excel) software.

If anybody has any suggestions I'm open to them all!


texts: How to convert an image file of lines into a lines shapefile in QGIS?

texts:
I need to convert the following image into a line .shp file. I can make the image into whatever type of image file it needs to be before converting. For context, I had started mapping an area by drawing over the map in photoshop but then came across QGIS and finished the project by mapping over the image using lines shapefile layers. I'm now trying to run the lines through a matlab code and the code can only use shapefiles.

I haven't used QGIS extensively.


texts: Arcobjects: how do I find if a line intersects a polygon?

texts:
I have one closed Polyline: "rectangle" and an open Polyline: "line".
Both Polylines may be Z aware.
What is the best way to find if the "line" intersects the "rectangle" in the XY plane?
That is I want a 2D (XY) test for intersection. I do not need the actual coordinates of any intersection(s).


texts:
I have MapServer setup to display data (both WMS and WFS), however I'm having a problem adding a default filter.

There's is an area on the map that users should not be able to retrieve features from. I have the coordinates of this area so basically, I need to setup a default spatial filter to exclude features that are contained within this area. Removing the data itself isn't an option.

The layers are all OGR layers.

I can construct a client side filter to do this, however this should really be a server-side setting.

How should I go about doing this?

Ignoring the spatial requirement for a second, I'm having problems getting even a basic filter to work. If I append the following to the URL of a WFS request, this filter works as expected:

(Note this filter has no actual relevance, it's just a simple filter so I can check if it's being applied)

&FILTER="<PropertyIsGreaterThan><PropertyName>ID</PropertyName><Literal>916</Literal></PropertyIsGreaterThan>"

However, if in the map file I set the following:

LAYER
   ...
   FILTER "<PropertyIsGreaterThan><PropertyName>ID</PropertyName><Literal>916</Literal></PropertyIsGreaterThan>"
   ...
END

Then I get no results from any query on that layer. Wrapping the filter in <Filter> tags has the same result.

I've also tried similar with the wfs_filter property, which doesn't seem to have any affect at all. (If I understand correctly, this is only used when using MapServer as a WFS Client?)

I'd appreciate any help.

Thanks :)


texts: How to create an atlas coverage layer from multiple lyers with multiple geometries in QGIS?

texts:
I have multiple layers (multi-point, multipart polygon, multipart lines). I want to create a hidden coverage layer based on all those to create an atlas of my QGIS project.

Obviously, I can't easily combine different geometries in one layer.

I ended up with the following workflow:

convert polygons to lines
convert lines to points
convert multi-points to single points
merge all single-point layers
create multipoint layer based on shared attribute (i.e., object
number)
use as hidden coverage layer

Since I will have to do that again in the near future, and this was a hell lot of manual work (and many layers): how would you do that properly?

Disclaimer: I read several Q/As which are somewhat related, but nothing fits, really. Most seem to deal with displaying multiple layers, not combining multiple layers in one coverage layer. Also, not being able to paste multi-part geometries in a single-part layer seems to be a currently persistent bug specific to Windows.


texts:
I could compile it with the following command on Macbook mavericks.

./configure --with-sfcgal=/usr/local/Cellar/sfcgal/1.0.4/bin

Then, I receive following results along with lists:

  PostGIS is now configured for x86_64-apple-darwin13.3.0

make gives an error:

  Can't fetch local revision (neither .svn nor .git found)
  Not updating existing rev file at 12547
  for s in liblwgeom libpgcommon postgis regress raster topology loader utils extentions; do \
  ...
  "sfcgal_triangulated_surface_triangle_n", referenced from:
  _SFCGAL2LWGEOM in liblwgeoma(lwgeom_sfcgal.o)
  "sfcgal_version", referenced from:
  _lwgeom_sfcgal_version in liblwgeom.a(lwgeom_sfcgal.o)
  (may be you meant: _lwgeom_sfcgal_version, _pg_finfo_postgis_sfcgal_version, _postgis_sfcgal_version )
  ld: symbold(s) not found for architecture x86_64
  clang: error: linker command failed with exit code 1 (use -v to see invocation)
  make[1]: * [postgis2.1.so] Error 1
  make: * [all] Error 1

Postgresql was installed using brew. I uninstalled Postgis brew installation in order to have SFCGAL functions compiled with Postgis.


texts: How to display table content using Qtableview after selecting a vector layer in Pyqt4?

texts:
I'm developing a QGIS plugin using QtDesigner and pyqt4.
I have a comboBox of my vector layers and I want to show the table of the selected layer.
I used this code for my comboBox:

def run(self):
    """Run method that performs all the real work"""

    allLayers = self.iface.legendInterface().layers()
    allLyrNames = [lyr.name() for lyr in allLayers]
    self.dlg.updateCombos(allLyrNames)

    # show the dialog
    self.dlg.show()
    # Run the dialog event loop
    result = self.dlg.exec_()
    # See if OK was pressed
    if result:
        # Do something useful here - delete the line containing pass and
        # substitute with your code.

        grdLayerName = self.dlg.getGridLayer()

        pass

And these two functions:

def getGridLayer(self):
    return(unicode(self.GridLayerCombo.currentText()))
def updateCombos(self, items):
    if len(items) > 0:
        self.GridLayerCombo.clear()

        for item in items:
            self.GridLayerCombo.addItem(item)


texts: How to Automatically Update ArcGIS.com Map Layer Data

texts:
We'd like to refresh layers in ArcGIS Online every five minutes with location data modifications occurring in our externally hosted SQL database.  We'd like to avoid purchasing ArcGIS for Server.

Any pointers for using REST to make the updates? Can we update multiple locations simultaneously? Can we push a JSON feed containing changes and deletions to ArcGIS Online every five minutes?


texts: How does Zonal statistics handle pixels that are on border of polygon?

texts:
I am using the plugin zonal statistics to calculate the sum of the raster cells within  polygons. How does this plugin handle pixels that are on the border of the polygon, are those cell values included in the statistics?

I am using QGIS 1.8.0.


texts: Google Earth Engine: Apply scaling function to band values of image

texts:
Situation
I would like to use a support vector regression model to predict a continuous outcome (percentage bareground cover). As predictors, I am using a number of different Landsat composite bands.
Here is a selection of band values:

As you can see, the band values are on very different scales, which leads me to believe that they should be scaled for the regression model.
Here is an example of a Landsat image with band values of very different scale (e.g. Band 9 vs. Band 10):
https://code.earthengine.google.com/?scriptPath=users%2Fjwiethase%2FTZ_restoration%3AUntitledFile2
Tried so far
It seems that a sensible scaling formula would be something like this:
scaled value = (x - mean(x)) / sd(x)

I know that I can apply some calculations to all bands of the image by using functions like .divide(), and calculate the average of a band using the reducer ee.Reducer.mean().
Question
How do I put it all together, to apply the scaling function to every value in every band of an image?
The desired outcome is an image with the same number of bands and band values, but each value scaled based on the function using the respective band average and standard deviation.


texts: CORS problem when trying to call WMS getFeatureInfoUrl using OpenLayers

texts:
I'm using the newest stable version of OpenLayers which is 6.3.1.

I have a problem when trying to call http request on an url genereted by the wmsSource.getFeatureInfoUrl function. It returns me an URL and when I execute this url directly on the newly opened tab on my browser I get the detailed result about the cliecked point on a map. But when my Angular app is calling this url I have a CORS error:

Access to XMLHttpRequest at 'http://q1.geopartner.ad/cgi-bin/qgis_mapserv.fcgi?map=/qgisdata/projekty/Q_postD.qgz&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetFeatureInfo&FORMAT=image%2Fpng&TRANSPARENT=true&QUERY_LAYERS=dzialki&LAYERS=dzialki&INFO_FORMAT=text%2Fhtml&I=126&J=31&WIDTH=256&HEIGHT=256&CRS=EPSG%3A3857&STYLES=&BBOX=2035059.4410645328%2C6966165.009797823%2C2074195.1995465432%2C7005300.768279833' from origin 'http://localhost:4200' has been blocked by CORS policy: No 'Access-Control-Allow-Origin' header is present on the requested resource.

My url won't be available for you cause it's only available behind the VPN connection.
Here is my code:

const wmsSource = new TileWMS({
  url: 'http://q1.geopartner.ad/cgi-bin/qgis_mapserv.fcgi?map=/qgisdata/projekty/Q_postD.qgz',
  params: { LAYERS: 'dzialki' }
});

const wmsLayer = new TileLayer({
  source: wmsSource
});

const view = new View({
  center: [ 2071206.0066185605, 6997895.536649533 ],
  zoom: 10,
});

const map = new Map({
  layers: [wmsLayer],
  target: 'map',
  view
});

map.on('singleclick', function(evt) {
  document.getElementById('info').innerHTML = '';

  const viewResolution = view.getResolution();
  const viewProjection = view.getProjection();

  console.log('viewProjection', viewProjection);

  const url = wmsSource.getFeatureInfoUrl(
    evt.coordinate,
    viewResolution,
    viewProjection,
    { INFO_FORMAT: 'text/html' });

  if (url) {
    console.log('url', url);

    const Http = new XMLHttpRequest();
    Http.open('GET', url);
    Http.send();

    Http.onreadystatechange = (e) => {
      console.log(Http.responseText);
    };
  }
});

The map is loading correctly, I can zoom in/out etc. And when I click on the map I have the problem which I described. How can I fix it?


texts: Cataloging small hardcopy map collections

texts:
We have a small collection of hardcopy maps, in the low hundreds, to catalogue and file into long term storage and are looking for lightweight guidelines and tools to do this in a way that will benefit future staff who aren’t going to remember “Oh yeah, the 1990s land negotiations! Those maps are in the cabinet behind the paper storage in the basement”.

We’re under strong pressure to “get those old pieces of paper out of here” to make way for desks and people and don’t have the time or budget to invest in a real library system. The solution of the day is to create a spreadsheet with a half dozen home-baked field names and a hyperlink to a photo or scan in the file-system. This would be a step or two above musty cabinets on the other side of the building but is a fragile solution.

Can you suggest something that might bridge the gap between home made spreadsheets and real catalogue created by a map librarian?


texts: How to style a LineString feature to be red colored?

texts:
I try to add Features with LineString geometries in different colors to a ol.source.Vector.

I tried to do this with the following code:

var myFeature = new ol.Feature({
   geometry : new ol.geom.LineString(coordinates),
   style : new ol.style.Style({
       stroke : new ol.style.Stroke({color : "red"})})
});
mySource.addFeature(myFeature);

The coordinates are correctly transformed. The resulting LineString Feature is displayed at the correct position in the map – but not in red.

What am I doing wrong?


texts:
I'm trying to implement a search box that uses my CartoDb data? My current code does not work.

Html:

<link rel="stylesheet" href="http://libs.cartocdn.com/cartodb.js/v3/3.15/themes/css/cartodb.css" />

My Map

 <form>
   <input type="Text" name="userInput" id="select-place" class="select-place">
   <input type="submit" value="Submit" onclick="processFormData();" />
 </form>

Js:

function processFormData() {
  var input = document.getElementById("select-place").value;
  return input;
}

function main() {
  cartodb.createVis('map', 'https://testuser.cartodb.com/api/v2/viz/123456789/viz.json', {
   title: false,
   description: false,
   tiles_loader: false
  })
  .done(function(vis, layers) {

  var userInput = processFormData();
  var sql = new cartodb.SQL({ user: 'testuser' });
    sql.getBounds("select * from testdb where name like ='%userInput%'").done(function(bounds) {
   map.fitBounds(bounds)
   });

var subLayer = layers[1].getSubLayer(0);
createSelector(subLayer);
 })
 .error(function(err) {
  console.log(err);
 })
 }
 window.onload = main;


texts: Can I delete the Styles_idx.dat file? in appdata/roaming?

texts:
Can I delete the styles_idx.dat file (in appdata/roaming - windows computer) - it is 125MB and causing trouble for my profile at the university, keeps saying my profile is full.

I have found the file, but did not want to delete it, if it would destroy my project?

Can I just shut the ArcMap program down, then delete the file and start up again?

What does this file actually represent?


texts:
I want to harvest from a WMS service like https://nsdiindia.gov.in/nsdi-portal/index.jsp in my geonetwork node and for which I am providing the details as given in the documentation

but I am getting the error like

  Unable to add a new entry
  Error : error
   Message : String index out of range: -1
    Class : StringIndexOutOfBoundsException

Can anyone help me with this?


texts: How to create a network with all the track points

texts:
I have created a sample network through the use of gps visualizer. Than I used the ogr2ogr function to import this data in my database.

Now I have a table which holds all the track points of these tracks.

I would like to build a network between all these points so that at the end I can calculate the shortest path from a point to any point. In all I have 28 track points where 26 of these points are distinct.

I would like some suggestions how can I determine the source and the target of each point in order to create a network.

Any suggestions please?


texts: Combining multiple DEMs into one file geodatabase

texts: I have multiple file geodatabases (approximately 240) each contain a DEM.  My issue is taking the DEM from each and putting it into a single file geodatabase without alot of mind numbing work.  Additionally the DEMs from each geodatabase is labeled DEM1.


texts: Using st_intersects caused mismatch when a marking (tagging) shape marked a polygon overlapping only on the edge

texts:
I have accidentally marked a landshape as being tagged "CU GO-M" because it seems to share an edge with a tagshape for "CU GO-M" but nothing else. I used the st_intersects method to join these shapes erroneously. Which method should I use? This might be a dumb question because I think it is st_within, but I'm worried that might kick out some false negatives (rather than a false-positive, ie this case).

To summarize, I think any interior overlap match should be marked as tagged. This case marks the tagging for only an exterior (edge) intersect.


texts: Debugging No Free Units (LUNGET) from RasterToASCII_conversion() in ArcPy?

texts:
I am converting approximately 20 years of satellite data to ASCII format (.asc) using ArcPy and rather than threading it, I chose to run four instances of the script in four different cmd windows. The code (slightly simplified) is as follows:

import arcpy, os
arcpy.env.overwriteOutput = True

years = range(1991, 2010, 4)
for year in years:
    for img in os.listdir("E:/tifs/{}".format(year)):
        if img[-3:] == "tif":
            output = "{}asc".format(img[:-3])
            print output
            arcpy.RasterToASCII_conversion("E:/tifs/{}/{}".format(year, img), "E:/ascii/{}/{}".format(year, output)

I have four of these running, with the start of the years range increased by 1 for each script copy. After going through about a year of the data, the script stops with a cryptic message of No free units (LUNGET), without any stacktrace:

This message does not turn up anywhere on ESRI help forums, stackoverflow, python forums, etc. The code runs fine without multiple instances of it running.


texts: Where can I find Road shapefile/e00/etc for the entire world?

texts:
I am hoping that there is a single source for this data, ie. 1 geodatabase. If not, can someone point me towards the sources that would make up the entire world's road transportation network?

I have come across references to the ESRI Digital Chart of the World, but I can't find any place to download it.

Ideas?


texts:
I am trying to iterate through features and perform a ZonalHistogram using each row. I am iterating using a SearchCursor. ZonalHistogram requires that I specify the zone field. When specifying the zone field of 'statsZone' I get the error below. I've ensured that this field does in fact exist.

import arcpy
import datetime
import os
import os.path
from arcpy import env
from arcpy.sa import *

newlayer = arcpy.mapping.Layer('D:/myshapefile.shp')
rasterLayer = "D:/myRaster.tif"
outputDirectory = "D:/myOutputLocation/"

with arcpy.da.SearchCursor(newlayer,["SHAPE@",'statsZone']) as cursor:
    for row in cursor:
        print(row[1])
        arcpy.gp.ZonalHistogram_sa(row[0], 'statsZone', rasterLayer, outputDirectory+row[1]+".dbf", "")


texts: Importing ECEF WGS84 into ArcGIS Desktop?

texts:
I have a file in ECEF WGS84 in meters.

I want to import these points in ArcMap 10.5 as point layer.

The file to import is in Excel format.
I can import them with Excel to layer tool. I can see all the points, but I can't get them on the map in the right position.

Here is the file header and first line:

file_name,file_type,sequence_id,timestamp,GPS_time(s),ecef(x),ecef(y),ecef(z),attitude(x)=roll,attitude(y)=pitch,attitude(z)=pan,frame_id
F:/DG_Paderborn-Borchen - Etteln_170620_Run1/2017-06-20_10-35-25/ladybug_panoramic_000000.jpg,JPG,0,28583339722,1181983021.1082721,3922148.9627709221,604199.61446186027,4976781.3346982058,-177.57368922265997,7.1055832760669553,-79.315623440207304,1068

Does somebody know a way to accomplish this?


texts: Example shapefiles containing PointM, MultiPointM, PolygonM or PolylineM data

texts:
I have written a Delphi library for reading shapefiles. Now I'd like to test it with some real data, but cannot find any shapefiles containing the "M" variety of data:

PointM, MultiPointM, PolygonM or PolylineM
and also "Z" data (PointZ, MultiPointZ, PolygonZ or PolylineZ) including the optional "M" data.

Any hints where I can get some test data of this kind?

https://gis.utah.gov/data/sgid-transportation/roads-system/ apparently links to some files but the ftp server times out on me.


texts: Select and load polygons intersecting another with PyQGIS (polygons stored in PostGIS)

texts:
I'm still a newbie at Python and I'm having issues selecting and loading polygons intersecting others with PyQGIS.
I managed to use PyQGIS to achieve that operation with shapefiles stored in my laptop but I haven't found a solution to proceed the same treatment on shapefiles stored in PostGIS.
This is what I came to so far :
from qgis.core import *
from PyQt4.QtCore import *
from PyQt4.QtCore import QSettings
from PyQt4.QtCore import QSettings
from qgis.core import QgsVectorLayer, QgsDataSourceURI
import processing

uri = QgsDataSourceURI()
uri.setConnection("localhost", "5432", "name of my database", "my user name", "my password")
uri.setDataSource("database schema", "table name", "geom")
couche = "PM"
layer1 = QgsVectorLayer(uri.uri(), couche, "postgres")
if not layer1.isValid:
    print ('Chargement de la couche: echec !')
else:
    QgsMapLayerRegistry.instance().addMapLayer(layer1)

uri.setDataSource("database schema", "table name2", "geom")
couche = "batiments"
layer2 = QgsVectorLayer(uri.uri(), couche, "postgres" )

subset_layer2 = processing.runalg('qgis:extractbylocation', layer2,layer1,['intersects'],0.0,None)

uri.setDataSource("rbal", "parcelles", "geom")
couche = "parcelles"
layer3 = QgsVectorLayer(uri.uri(), couche, "postgres" )

subset_layer3 = processing.runalg('qgis:extractbylocation', layer3,layer1,['intersects'],0.0,None)

When I run the code I don't get any error but the treatment is loading forever without returning any results...

EDIT:
With that script I've managed to select all the polygons contained in layer2 that intersect layer1 but I do have to load layer2 to do so.
Is there a way to select and load those intersecting polygons without loading the whole layer in the first place ?
uri = QgsDataSourceURI()
uri.setConnection("localhost", "5432", "name of my database", "my user name", "my password")
uri.setDataSource("database schema", "table name", "geom")
couche = "PM"
layer1 = QgsVectorLayer(uri.uri(), couche, "postgres")
if not layer1.isValid:
    print ('Chargement de la couche: echec !')
else:
    QgsMapLayerRegistry.instance().addMapLayer(layer1)

uri.setDataSource("database schema", "table name2", "geom")
couche = "batiments"
layer2 = QgsVectorLayer(uri.uri(), couche, "postgres" )
if not layer2.isValid:
    print ('Chargement de la couche: echec !')
else:
    QgsMapLayerRegistry.instance().addMapLayer(layer2)

feats_layer1 = layer1.getFeatures()
for feat in feats_layer1:
    geom_layer1 = feat.geometry()
feats_layer2 = layer2.getFeatures()

for feature in feats_layer2:
    if feature.geometry().intersects(geom_layer1):
        layer2.select(feature.id())


texts: Generate N distinct colors which (alpha) blend in a nice manner when overlaying an orthophoto or a map

texts:
I have N polygons (let N = 3 for making it simple) to display onto a base map in an OpenLayers application (but this can also naturally be applied to other GIS softwares as well).
For the moment, I manually picked the color for each polygon that suits my eyes, adding an alpha value to the fill so that I can keep the details of the underlying base map as much as I can (but this obviously is a trade-off between knowing which polygon overlays the map and what are the terrain details under).
Here is the code that generates the colors:
let alpha = 0.9;
const colors = [
    [136, 0, 255, alpha],
    [255, 243, 138, alpha],
    [143, 251, 255, alpha],
    [255, 161, 179, alpha],
]

let borderColors = [];
let fillColors = [];

for (let row of colors) {
    let rowstr = JSON.stringify(row);
    // weired and probably useless way of converting each row of the colors array:
    borderColors.push(rowstr.replace(/\[/g, 'rgba(')));
    fillColors.push(rowstr.replace(/\[/g, 'rgba(').replace(/0.9]/g, '0.2)'));
}

Which I then use in the style definition of a vector layer, in a loop over my N polygons:
style: new ol.style.Style({
    stroke: new ol.style.Stroke({
        color: borderColors[i],
        width: 2
    }),
    fill: new ol.style.Fill({
        color: fillColors[i],
    }),
}),

Results are as follow on an orthophoto (which is more critical than a map because of its main 'greenish hues'):

On the left, the original colors of the base layer. Then, going toward the right, a cyan polygon, which then has a yellow one on top, and finally a pink-red one over the two previous ones.
Therefore, each time one or more polygons overlay an original pixel p, this pixel is getting moved in the colors space by some amount. I would like to be able to distribute, on average for a standard orthophoto (I have the feeling that it wouldn't be useful to do some calculations, which would also be very time consuming...), every possible shift in the color space so that all combinations of overlap between polygons

are as discernible as possible from each other,
keep the underlying imagery contrast, and
pleasing to the eyes.

I have seen Generating random muted/pastel colours using QGIS  which was asked recently, and is great and inspiring, but in my case, the main difference is that the N polygons have (2^N)-1 different overlapping possibilities (i.e. 7 my case with 3 polygons). So I need to pick 3 colors which 7 combinations can be mathematically described to render a nice cumulative stack (using the alpha channel) of colors on the screen.
Maybe part of the solution could come from the classic overlay modes offered by any imaging software, like in Gimp where you can choose a blending mode:

In this situation I wonder if that would be possible in OpenLayers / JS in general?
Which blending mode is used by default on OL, and which mode would give better results?
Another idea would be to create a number of individual layers for each of the different overlapping possibilities (7 in this case), and fall back on the related question. But I would like to avoid this because each polygon has intrinsic attributes that I must keep attached to it.


texts: Stretching of geometry when rasterizing vector data on Google Maps

texts:
I'll start with a quick backstory on what I am trying to achieve.  I'm working on an offline Android GIS tool and people are throwing giant datasets at it.  It converts these into Google features (V2) and in the case of shapefiles or KML's with tens of thousands of features it overloads the rendering capacity and memory, making the map choppy.

I had the idea of rasterizing these datasets on the fly.  My approach is a little naive in that I am simply converting the Lat/Longs to image space, drawing the features on a canvas and then throwing the Bitmap on the map as a GroundOverlay.  Over small areas it is pixel perfect:

However over a much larger geographical area there is a pronounced vertical warping of the geometry:

The boundaries of the rasterized image line up perfectly with the datasets but the internal features warp vertically away from the equator.  I'm guessing there is some sort of projection issue at play but I am lost.

Could someone explain what I am seeing and how to correct the alignment?


texts: Determine the flow direction from a point in a DEM

texts:
I had a quick look of the other questions about hydrology using a DEM (SAGA/ArGIS) and I guess I need to calculate something different from a traditional flow accumulation/channel network/catchment basins for an hydrological model.

I'm working with pipe network systems and when the capacity of system is over, the water is coming out from a single or more manholes. What about the water, where is going?

The manhole is a single point in a DEM, is there an easy tool to try to assess the most probable flow path of the water over the DEM?
(Without calculate the stream network for the whole DEM until accidentally a stream cross the manhole)

I'm mainly using SAGA and ArcGIS. I read and tried most of the tools in Saga (SAGA - Tools - Hydrology) but none seems to be what I'm looking for. They are all trying to calculate a channel network for the whole DEM when I need just a line from a point as output.
The only one useful is the Downslope Area but it is interactive (I have to select the starting point with the mouse). The not interactive version seems to be released in a new version of SAGA (Feature Request -SAGA).

My questions are:

Are there similar tools in SAGA?
Is there a way to start the Downslope Area from many different points?
Is there a similar tool in ArcGIS?


texts: getParameterInfo datatype for creating new feature class

texts: I'm having trouble finding a suitable datatype that allows the user to browse to a .gdb and create a name for a new feature class. DEFeatureClass wants an existing FC to write to. Using CreateFeatureclass_management to do so before the user browses to the .gdb and enters the name doesn't seem to work. This is for a Python Toolbox in 10.2.2


texts: Iterating files in ModelBuilder of ArcGIS Desktop 9.3?

texts:
I want to use around 100 files as input in a model I created in ArcMap 9.3

Do to so I followed the esri help page (http://webhelp.esri.com/arcgisdesktop/9.3/index.cfm?TopicName=An_overview_of_model_iteration) and right-clicked my input variable to switch it from "A single value" to "A list of values". When I run my model from the toolbox it only offers me list where I have to input each of my files separately. It is not possible to select all my 100 files at once like I can do in ArcMap 10.x when using the "Iterate Feature Classes" tool.

Is there a way to select a bunch of files as input for a model in ArcMap 9.3 without selecting each file separately?


texts: I am using Mosaic to New Raster for AvgMin temp for each month, I need to reference from which month any given pixel in the output was populated from

texts:
In ArcMap, I am working with climate data, AvgMin for each month. Using Mosaic to New Raster, and Mosaic Operator set to minimum I can isolate the absolute lowest value over the 12 input rasters, however, I need some way to record which month the value was populated from.

Ideally, I would like to see one column with a temperature value, and another with a value ranging from 1-12 to reflect which month the temperature value originated from.


texts: QGIS 3.2 - How can I scale polygon about his center?

texts: How can I scale a polygon about his center? I find only transform function, but the scale function here can only move the polygon with the scale factor and don't enlarge it.


texts: Problem dry fields- urban areas on supervised classification in QGIS, Landsat 7

texts: I am trying to make a supervised classification on a Landsat 7 image using the QGIS 3.16.1 version. I am using the SCP plugin to do then classification.  I am doing the classification using the RGB band 4-3-2 and the maximum likelihood classification.The problem I have is that in the final product the bare or dry fields appear like urban areas. Does anyone know how to solve this?


texts:
I am experimenting drawing curved line from Origin to Destination in Cartodb.
I was able to do so once, however, failed the second and third time.

Could anyone offer me some insights into the following code and help me debug?

I am using CartoDB and start with a table of centroid points for all the countries.

Step 1:CHANGE THE REFERENCE SYSTEM
I used the following code

INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) values ( 953027, 'esri', 53027, '+proj=eqdc +lat_0=0 +lon_0=0 +lat_1=60 +lat_2=60 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs ', 'PROJCS["Sphere_Equidistant_Conic",GEOGCS["GCS_Sphere",DATUM["Not_specified_based_on_Authalic_Sphere",SPHEROID["Sphere",6371000,0]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Equidistant_Conic"],PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",0],PARAMETER["Standard_Parallel_1",60],PARAMETER["Standard_Parallel_2",60],PARAMETER["Latitude_Of_Origin",0],UNIT["Meter",1],AUTHORITY["EPSG","53027"]]');

Step 2:CREATE THE LINES

update world_borders_coarse_2013
set the_geom=
(select st_transform(
  st_segmentize(
    st_makeline(
      st_transform(world_borders_coarse_2013.the_geom,953027),
      st_transform(b.the_geom,953027)),100),4326)
 from world_borders_coarse_2013 b
 where b.name='China'
 )

It should give me the curved lines between China to all the other places. However the reality is that I have the "the_geom" column as "line", but nothing shows in the "map view".

Would anyone give me some suggestion to change the code?


texts: I took a tif file and cropped it to a portion using Raster>Extraxtion>Clip Raster by Extent. My project name is OBIA. After cropping the image is saved to this project. How do I download this image?


texts: No maps visible when downloading something with CartoPackageManager

texts:
I am facing such problem: while some package is downloading with Carto PackageManager, I can’t see neither other downloaded offline maps nor online maps of any other place.

In other words, if there is some package that is pending or downloading, I can't use any map,
I can see only gray background from offline or online vector sources in any location, but as soon as map packages are fully downloaded the map appears in a regular way.

In debugger in both cases no differences, I easily pass over lines (this is full code of mapView initialization) :

        MapView = new MapView();
        ConstraintUtils.StretchView(MapContainerView, MapView);

        var baseLayer = new CartoOnlineVectorTileLayer(CartoBaseMapStyle.CartoBasemapStyleDefault);
        baseLayer.Language = MapElementsConstants.LocalizationLanguage;
        MapView.Layers.Add(baseLayer);
        _source = new LocalVectorDataSource(MapView.Options.BaseProjection);
        _layer = new VectorLayer(_source);

        MapView.Layers.Add(_layer);
        MapView.SetZoom(zoom: MapElementsConstants.MapControlHighZoomingCoeff, durationSeconds: 0);
        MapView.UserInteractionEnabled = false;

Is this feature of CartoSDK or it may be a bug in my application?

All iOS versions (11.2 in particular), Carto NUGET Package version : 4.0.2


texts: What does n mean in this formula: n = pow(2, zoom)? I using this when I calculate tiles.


texts: Changing info div based on baselayer selection

texts:
I would like to change the map logo in a info div based on the baselayer selection.
From what I understood, I have to use events but I do not know how.
Here is my code:

var info = L.control({position: 'bottomleft'});

info.onAdd = function (map) {
    this._div = L.DomUtil.create('div', 'info'); // create a div with a class "info"
    this.update();
    return this._div;
};

map.on('baselayerchange', function(eventLayer) {
if (eventLayer.name === 'MapboxDark') {
    info.update = function (props) {
            this._div.innerHTML = '';
            };
}
if (eventLayer.name === 'MapboxOutdoor') {
    info.update = function (props) {
            this._div.innerHTML = '';
            };
            info.addTo(map);
}

});
info.addTo(map);