texts: How to read the relevant join information out of a MXD file (above all the name of the destination table attribute)

texts:
Goal: Read all joins within a MXD file and its relevant information. This should include: origin layer name, origin attribute name, destination table name, destination table attribute.

To solve this, I used the code snippet from the following post:
Seeing all existing joins in any MXD layer in one place?

It all works so far. I am able to read the origin layer name and attribute name. However I'm only able to read the table name of the join table (destination table name), not the name of the attribute which is used to build the join (RelationshipClass.DestinationForeignKey is always empty). The joins are between FeatureClasses and geodatabase tables (enterprise and file-based). I'm probably missing something...

My code looks like this:

# Main Function
def get_lyr_join_fields(self, display_table, layer_name):
    return_list = []
    while self.__try_cast(display_table, db.IRelQueryTable):
        source_table_name = layer_name
        rel_query_table = self.__try_cast(display_table, db.IRelQueryTable)
        rel_class = rel_query_table.RelationshipClass

        # Join Source Tabelle
        source_table = rel_query_table.SourceTable
        source_dataset = self.__try_cast(source_table, db.IDataset)
        source_class = rel_class.DestinationClass
        source_att = rel_class.OriginForeignKey

        # Join Destination Tabelle
        dest_table = rel_query_table.DestinationTable
        dest_class = rel_class.OriginClass
        dest_table_name = dest_dataset.Name
        dest_att = rel_class.DestinationForeignKey # <--- ALWAYS Empty!

        # Build retrun-List
        return_list.append("{}:{},{}:{}".format(source_table_name, source_att, dest_table_name, dest_att))

        # Get other joins
        display_table = source_table
    return return_list

# Cast-Method (for COM-Objects)
def __try_cast(self, obj, interface):
    try:
        new_obj = obj.QueryInterface(interface)
        return new_obj
    except:
        return None

While debugging the code I get the following output for the RelationshipClass (join Attribute of the origin is read but not for the destination):


texts: How to save an existing LYR file to use a relative path *and* have the LYR be v9.3 compatible?

texts:
Not sure if this is possible - I am writing a python script to convert LYR files to use relative paths and save the LYR as a 9.3 LYR file.  I can do each one individually, but not together.

...

# save existing LYR to a new file and set to RELATIVE path - works
arcpy.SaveToLayerFile_management(lyr, relFile, "RELATIVE")

# saveACopy to v9.3 LYR file - works, but does not maintain RELATIVE path
lyr2 = arcpy.mapping.Layer(relFile)
lyr2.saveACopy(93File,'9.3')

...

It doesn't seem to matter if I reverse the order of these calls, I still get v9.3 LYRs with absolute paths and v10.x LYRs with relative paths.  I have run this in 10.0 and 10.1 and see the same results.


texts: Are interpolation methods the same in all software, produce the same result ( with same settings) or does it differ from software to software?

texts:
I was wondering if interpolation algorithm has been changed? produce the same result in all software so it is not software based. This would be if the settings are exactly the same of course.

Let's say that I use two different software to produce a nearest neighbor interpolation with the same data. Will it be the same result or is the algorithm a bit different for every software?


texts: Is there a simple way to detect what generation a MrSID file is?

texts:
I just found out there are a few types of MrSID files - generations MG2,MG3,MG4. As I have a user which uses Bentley Microstation V8 and apparently can see only some of the generations (and thus I need to convert between the different generations)- I wanted to know if there is a quick way to detect what generation a SID file is (I have AutoCAD 2007, ArcGIS 9.3.1)?

Thanks in advance.


texts: Free And Open Source Software Tools to orthorectify Pleiades imagery

texts: Are there Free And Open Source Software Tools to orthorectify Pleiades, Worldview And Quickbird imagery as we do with COST tools such as ERDAS Imagine, PCI, ENVI, etc.?


texts: Is there Open Source equivalent to ArcGIS Integrate (for topology correction)?

texts:
I need to transform a set of polygons so that it forms an exact spatial partition/tesselation (no gap - no overlap). What are the existing open-source solutions (software or library) to perform this correction? Ideally, I am looking for open-source versions of the ArcGIS integrate.

  Integrate is used to maintain the integrity of shared feature
  boundaries by making features coincident if they fall within the
  specified x,y tolerance. Features that fall within the specified x,y
  tolerance are considered identical or coincident.


texts: Polygon To Centerline produces weird end in ArcGIS Pro?

texts: When using the 'Polygon To Centerline' tool on a river polygon it creates a centerline with a very nice end and a very weird end (see pictures). Any ideas what the problem could be? I would like to have 2 nice ends. I use ArcGIS Pro 2.2.0 with an advanced license.


texts:
I just added the GeoJSON vector layer in Leaflet for my hybrid mobile app using PhoneGap, the problem is that when I add the GeoJSON file of size larger than 2MB the application crashes.

What is the problem can larger GeoJSON files not be loaded this way?

The code for loading GeoJSON file in my map window is below:

var roadStyle = {
     color: "#00ff00",
     weight: 3
};
var roads = [geojson]; //GeoJSON data is added in roads variable
L.geoJson(roads, {
 style: roadStyle
}).addTo(map);


texts: Selecting Features in ArcGIS based on their relation other layers

texts:
I am extremely new to GIS mapping and I'm trying to do what I believe would be relatively simple. I have two layers: one contains data on school locations in Chicago and the other contains data on child abductions; both layers are point shape. What I would like is to only show the schools that are within one mile of a child abduction in order to reduce the number of features shown on the map.

I apologize if this is a duplicated question, but the 'Similar Questions' box didn't suggest the type of question I'm looking to have answered.

Thanks in advance for your help.

Technical details: I'm operating on a Windows 7 OS using ArcEditor Single Use evaluation trial software. I'm somewhat familiar with Python (but not overly familiar) if needed.


texts: Create polygons from raster with only a selected amount of features (not each pixel a different polygon)?

texts:
I am trying to convert a raster into different polygons. The raster is based on a DEM and only contains three classes (each pixel has value 1, 2 or 3, obtained by using the raster calculator).

I want to create three polygons out of this (so just one shapefile with three features) but when I use the polygonize processing tool the result is a different ID for each pixel and as an attribute the value 1,2 or 3.

However in a later step I want to assign an average point value to each polygon and in this way it is not working.


texts: CSV attribute table to polyline shapefile

texts:
I have a .csv file that is an exported attribute table that I would like to reconstruct as a polyline shapefile (for reasons beyond my control I cannot access the original shapefile from which the .csv attribute table was exported, nor do I have access to metadata). To reconstruct the polylines, I must:

Import the .csv file,
Display the point data, then
Connect the points in the correct order to make lines.

I am able to complete  steps 1) and 2), but am having trouble with step 3).

Here is what the exported .csv attribute table looks like (hiding a few irrelevant columns):

It seems like each row represents a line and that the field "the_geom" contains the lat/long of the points that make up the line. I've tried unsuccessfully to use the xy to line and points to line tools to reconstruct.

How can I reconstruct these lines?


texts: Receiving submodel or tool name as input using QGIS or ArcGIS ModelBuilder?

texts:
I would like the user to select the tool (or submodel) to be used from among a list of names.
Can I do this in QGIS or ArcGIS ModelBuilder without having to write a python script?

Actually, what I would like is the user to input in a table as below (like in the "union" tool in ArcGIS)

First column: a raster layer. Second column: the submodel or tool selected from a drop down list. Then a couple of columns for inputting values to be passed on to parameters in the selected submodel. I would also like the user to be able to add as many rows (i.e. rasters) as necessary.


texts: Does PostGIS have separate implementation of GiST index from PostgreSQL?

texts: I want to know if PostGIS has a separate implementation of GIST index or it reuses the GIST index implementation provided by Postgresql?


texts: how to transform DB stored EPSG data to Geography data type in sql server

texts:
I have a 3 GB database in SQL Server, with certain points stored in float values. The points are stored in EPSG 3765. I would prefer to have them saved in the Geography datatype. I can not store them directly, because the points are stored like this 1734994.3283045683    5740529.6114465743 ( x, y )

Geography requires values between 90 and -90 as well as 180 and -180.

How can I do this? I've been Googling for the last 4 days and can't find a solution. I wouldn't mind writing any code (C#, Python, SQL, etc)

I have no shapefile or anything else, just the data in the DB.


texts: Making summary table of continent wise contribution of research papers in literature using QGIS

texts:
I want to generate a summary table. This table will consists of contribution of each continent in research literature.
I have a CSV file that consists of latitude and longitude information of each affiliation. I am attaching link of a csv file.
I want that kind of table. How can I get it?


texts:
I am a newbie to R and I want to divide a country map into grids so that I can come to a decision according to the latitude and longitude values which grid am I residing.

I am looking for an example (R code) how to this.


texts: Getting p-values from mgcv GAM result and calculating their average

texts:
I'm doing many GAM model calculations (mgcv). Their output looks like the attached image. What I want to do is 1) to get only p-values from the output and 2) calculate their average. For example, in the case of the image, I want to get 0.0477 and 0.1518 and get their average (0.09975). I can do one calculation manually but I have many models and want to do it automatically.

I did look at pages like below but I get p-values different from the summary.

(added)

The format of GAM model I used is:

model1 <- gam(dependentx  ~ s(independent1) + s(independent2), data=datatable, method = "REML", bs='cr', family=Gamma(link=log))

https://stackoverflow.com/questions/53389900/extract-p-value-from-gam-check-in-rhttps://stackoverflow.com/questions/53389900/extract-p-value-from-gam-check-in-r


texts:
For my bachelor thesis in biology I´m working with QGIS 3.14 (Pi) and tried to create a layer where the slope aspect/expositon/exposure is presented. I´ve got a DTM1 on hands. There is this function of Qgis: raster --> raster terrain analysis --> aspect  - with which I was able to create a nicely looking map. The function gives you an azimuth number (between 0-360), which I assigned to the corresponding compass direction. 0-22,5 = N; 67,5 = NE etc.
BUT what I get out of there is certainly not the aspect what it is in nature.
It kinda seems like north and south are inverted sometimes. East and west seem to be ok.
I attach some pictures to show you what I mean.

With only hillshade on you can see, that the slopes my program is pretending to be north in truth are NOT north. Rather south/east.

Do you have any Idea what could have gone wrong?
This is my GDALinfo:


texts: Changing nearest neighbor parameter in Dzetsaka plug-in

texts:
I am using the Dzetsaka classification tool plug-in in QGIS 3.10.6. I am testing out different methods such as Random Forest and KNN. I was wondering if there is a way to change the nearest neighbor parameter when using KNN? I do not see any options to do so.
Any suggestions?


texts: Please suggest me how can I fill the gap of SLC off Landsat 7 data band 8 that is panchromatic band. Also suggest me is there any method to fill the gap of all other bands too.


texts:
I have this script I am trying to use to search for a value within a CartoDB table. I can't seem to figure out how to reverse the order of the lat & long coordinates so they appear as lonlat... instead of latlon

Here is what I am trying to do:

function parcelSearch(layer) {
        var sql = new cartodb.SQL({ user: 'dchadney' });
        var $options = $('#searchParcel');
        $options.click(function(e) {
            var streetNam = document.getElementById('parcelid').value;
            if(streetNam === "") {
                alert("Please enter the parcel ID");
            } else {
                query = "SELECT * FROM parcels_sa WHERE parcel_id ILIKE '%" + streetNam + "%'";
                console.log(query);
                layer.getSubLayer(0).show();

                sql.getBounds(query).done(function(bounds) {
                    latlon = bounds[0];
                    console.log(latlon);

                    $('.latlongit').val(latlon);
                   // map.fitBounds(bounds);
                });

                    $('#step1').modal('hide');
                    $('#step1A').hide('fast');
                    setTimeout(function() {
                    $("#submitModal").modal({
                    backdrop: 'static'
                    });
                    }, 7000);
                    $('button#step3indicator').addClass('active');
                    $('.line2').css('width', '42%');
            }

            // change the query in the layer to update the map
            layer.setSQL(query);

        });

Seems to me that latlon is a function coming from leaflet? How do I modify this?


texts:
I'm new to working with gis software. I've always been a field guy who sends all his data in to the office for processing.
However I am trying it myself and have ran into problems pretty early on.

Basically what I want to do is take the raw data from my Garmin and convert it into a CSV which I can then import into QGIS for further processing. I was told that DNRGPS would be the program to use, but I can't seem to find any version of this or any other that will work with ios processing system.

I am open to any and all suggestions.


texts: Removing land based cells from sea surface temperature grid by using shoreline shapefile?

texts:
I am a second year Geo student.

I was tasked with creating a sea surface temperature grid from a colour profile, which I have completed, but I am required to drop the land based cells through the use of a shoreline shapefile.

I have dissolved the shoreline (which I am hoping is the right move), but can't seem to find a way to move forward.


texts: Using UpdateCursor to populate new field based on conditional statement not updating field

texts:
I am attempting to populate a new field with a "Y" or "N" based on whether two integer fields are equal.  My code seems like it should work, but when I run it, it does not update the field (the field is still ).  I have double checked that the inputs are correct and I have reloaded the table to make sure it was not an updating issue.

inTable = os.path.join(env.workspace, "pivotTableTEMP")
with arcpy.da.UpdateCursor(inTable, ["dmready", "total_records", "ALLREADY"]) as cursor:
    for row in cursor:
        if row[0] == row[1]:
            row[2] == "Y"
            cursor.updateRow(row)
        else:
            row[2] == "N"
            cursor.updateRow(row)

The dmready and total_records fields are both integers and the ALLREADY field is string.

Vince's advice worked perfectly for the above code.  However, I added another conditional statement and it is not filling properly.  Here is my new code:

inTable = os.path.join(env.workspace, "pivotTableTEMP")
with arcpy.da.UpdateCursor(inTable, ["dmready", "total_records", "dm_stat", "ALLREADY"]) as cursor:
    for row in cursor:
        if row[0] == row[1] and row[2] == "Ready for DM":
            row[3] = "Y"
            cursor.updateRow(row)
        else:
            row[3] = "N"
            cursor.updateRow(row) # could also delete here if

The result of this code is that it fills all records of the ALLREADY field with "N" when there are several that fit the conditions for "Y".  Any thoughts?


texts:
I converted ecw file into jpg file with IrfanView, to display the jpg as base layer in openlayers.
I understood that I need to create a jgw file to display the jpg file.

I found this example:

 4.23214625853148 - A
 0.00000000000000 - D
 0.00000000000000 - B
-4.23214625853148 - E
 3404018.70881921 - C
 5819863.55539414 - F

for F - i know the left up corner coordinates(in Israel TM projection), i need other projection?

I did not understand what is the other parameter:
A is the x scale length in meters? and E its the negative?
And about D and B I have no idea.

The file is little bit big: about 300 MB, its will be a problem in openlayers?


texts: How to find the extent of (many) EPSG CRSes?

texts:
I want to map the geographic extents / bounding boxes of a lot of EPSG CRSes. At least, those with a rectangular geographic extent.

For example, take epsg:2004, visiting that page shows the extent (both projected and geographic).

There are several sites with this info, which make it trivial to find the extent for a single CRS at a time

spatialreference.org
epsg.io, including via the QGIS "EPSG CRS" Plugin
epsg-registry, in the 'Area of Use' results.

Is there an easy way to extract this information off-line using QGIS, proj4 or other FOSS tools? Preferably without needing to parse the WKT/proj4 strings (which are easy enough to find, they're sitting in any postgis-enabled database)

Online queries (to an API or web service) would be good too.

Other things I've looked at:-

I don't see anything in the QgsCoordinateReferenceSystem class for this, so I assume that qgis doesn't do any bounds checking and throws all the points at proj4, which will complain (with a "forward projection failure") if any points fall outside the valid bounds.
Also, pyproj (the python wrapper for proj4) doesn't seem to be able to extract this info, it's really just for transforming between CRSes.

EDIT

Following BradHard's suggestion, I registered on epsg.org and built the postgres database using the scripts provided.

Then it was simple matter to join the CRS to the area...

SELECT
  epsg_area.area_name,
  epsg_area.area_of_use,
  epsg_area.area_south_bound_lat,
  epsg_area.area_north_bound_lat,
  epsg_area.area_west_bound_lon,
  epsg_area.area_east_bound_lon,
  epsg_area.area_polygon_file_ref,
  epsg_coordinatereferencesystem.coord_ref_sys_code,
  epsg_coordinatereferencesystem.coord_ref_sys_name
FROM
  public.epsg_area,
  public.epsg_coordinatereferencesystem
WHERE
  epsg_coordinatereferencesystem.area_of_use_code = epsg_area.area_code
  and epsg_coordinatereferencesystem.coord_ref_sys_code = 27700;


texts: Is it required to use BBOX strategy with WFS protocol in OpenLayers 2?

texts:
is it a must to use BBOX strategy with WFS protocol in Vector layers ?

I cannot display this layer without BBOX strategy:

vLayer = new OpenLayers.Layer.Vector("Test", {

    styleMap: new OpenLayers.StyleMap({
        "default": style
    }),
      maxExtent: new OpenLayers.Bounds(-180, -90, 180, 90).transform(geographic,mercator),
      sphericalMercator: true,
      protocol:  new OpenLayers.Protocol.WFS({
      url: 'http://myDomain/geoserver/wfs',
      featureType: "Layer_1",
      featureNS: "http://mapmap.org",
      srsName: "EPSG:3857",
      version: "1.1.0",
      extractAttributes: true,
      isBaseLayer: false,

        })
      });

I am not getting any request in firebug.
What am I missing here ?

EDIT:

URL request in firebug:

url=http://myDomain/geoserver/wfs

XML Post:

<wfs:GetFeature xmlns:wfs="http://www.opengis.net/wfs" service="WFS" version="1.1.0" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
<wfs:Query typeName="myWorkSpace:Layer_1" srsName="EPSG:3857" xmlns:myWorkSpace="http://mapmap.org"></wfs:Query>
</wfs:GetFeature>

and sometimes the response came with error:

<?xml version="1.0" encoding="UTF-8"?>
<ows:ExceptionReport version="1.0.0"
  xsi:schemaLocation="http://www.opengis.net/ows http://myDomain/geoserver/schemas/ows/1.0.0/owsExceptionReport.xsd"
  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:ows="http://www.opengis.net/ows">
  <ows:Exception exceptionCode="NoApplicableCode">
    <ows:ExceptionText>java.lang.IllegalStateException: A transformation exception occurred while reprojecting data on the fly
A transformation exception occurred while reprojecting data on the fly
Latitude 91°00.0'N is too close to a pole.</ows:ExceptionText>
  </ows:Exception>
</ows:ExceptionReport>


texts: Splitting a line layer using another point layer

texts:
How to split a line layer using another point layer in qgis?

I have one point layer and one line layer. I want to split the line layer between two points (those two points belongs to point layer) using qgis. how can I perform this splitting using qgis? I used v.split for splitting
the line layer, but I am not getting any output.


texts: Using os.walk and addDataFromPath to add multiple shapefiles from directory to map

texts:
I am attempting to add all the shapefiles from a single directory into an ArcGIS Pro map. Here is my code:
import arcpy
import os

path = # directory path
aprx = arcpy.mp.ArcGISProject("CURRENT")
map = aprx.listMaps()[0]

for dirpath, dirnames, filenames in os.walk(path):
            for file in filenames:
                if file.endswith("shp"):
                    map.addDataFromPath(os.path.join(dirpath,file))

But this only manages to add the very first shapefile that it finds, and then provides this error:
Traceback (most recent call last):
  File #script location, line 12, in <module>
    map.addDataFromPath(os.path.join(dirpath,file))
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\_mp.py", line 1771, in addDataFromPath
    return convertArcObjectToPythonObject(self._arc_object.addDataFromPath(*gp_fixargs((data_path,), True)))
RuntimeError
 Failed to execute (AddShapefiles).

However, when I replace map.addDataFromPath(os.path.join(dirpath,file)) with arcpy.AddMessage(os.path.join(dirpath,file)), the geoprocessing message does actually list out every shapefile within the directory.
What am I doing wrong?


texts: Mapnik - Weird tiles rendering when used with spatialite

texts:
I use Mapnik in an GeoDjango App with Spatialite 4.0 as backend. There is something weird with the tiles that I suspect might be related with a bad handling of mbr with spatialite.

If you look at those two images, the first is the tiles rendered in OpenLayers with a map.zoomToExtent on the features. If you look at the triangle in atlantic ocean you will see that the feature is not rendered entirely by Mapnik. The second image show the entire feature not tiled.

If I zoom in, the triangle continue to not be tiled correctly. If I zoom out a lot it finally appears.

I use this query in Mapnik:

query = '(select ' + geometryField + ', "id" from "shapefile_feature" where' + ' shapefile_id=' + str(shapefile.id) + ') as geom'
datasource = mapnik.SQLite(file='C:\mygeosite\sqlite3\sql3.db',
                                table=query,
                                srid=4326,
                                geometry_field=geometryField)

featureLayer = mapnik.Layer("featureLayer")
featureLayer.datasource = datasource


texts: Shapely affine transformations with bounds of GeoDataFrame

texts:
I understand that I can use Shapely's affinity.translate to translate a geometry like the below:

from shapely import affinity
from shapely.geometry import Polygon

poly = Polygon([(0, 0), (1, 0), (1, 1)])
poly2 = affinity.translate(poly, xoff=1.0, yoff=0.0)

How can I apply this transformation on the bounds of each polygon geometry in a GeoDataFrame and automatically update the position of the geometry as well?


texts: Easy way to see the requests QGIS makes to a remote server?

texts:
I'm attempting to debug an issue with a remote ArcGis Server's WMTS as seen from GeoServer. I can see the requests GeoTools is making (because they throw an error and get logged). However the service works in QGIS (3.0 if that makes a difference) and I would like to know what the difference in the requests is.

I know I can set up wireshark or some such but I hoped there was a way of getting QGIS to log the URLs being requested directly via a flag or setting?.


texts: Cannot Configure Web Adapter After ArcGIS Server Upgrade

texts: My company just upgraded from Server 10.3 to 10.5. We are trying to configure the web adapter and keep getting errors about invalid tokens after entering in our FQDN "http://METCONAPP.METCON.DOMAIN:6080" We got this domain name from the original web adapter config file. Our default web site adapter in IIS was called "arcgis" and the installer wouldn't let us name the new adapter "arcgis," because the folder was still in the webroot. We tried to delete this web adapter from the webroot folder and start over, but we still haven't made any progress. WE're not sure how to proceed. We need to get our website and ArcGIS server back online ASAP.


texts: How to edit and use text produced by Feature InfoTool?

texts:
How can I access and use text from Feature InfoTool dialog. I'm working with JMapFrame. I found this snippet:

         InfoTool tool = new InfoTool() {
         @Override
         public void onReporterUpdated() {
             String text = getTextReporterConnection().getText();
             System.out.println(text);
         }
     };

but I don't see how to link this to my JMapFrame. So basically I don't know how to get InfoTool from JMapFrame.


texts: ArcGIS Server Java Edition - Deploy App to Embedded Tomcat

texts:
I am more experienced with the .NET side of the house, so please excuse me if this is a simple question.

Can I deploy an non-ArcGIS/ESRI app to the embedded Tomcat instance shipped with ArcGIS Server 10 SP5?  I have created a simple JAX-RS Jersey REST web service that pulls lookup values from MSSQL db to feed various widgets.  As this is (currently) the only component outside of ArcGIS Server ecosphere - I am hoping that I can use the embedded tomcat instance vs installing a more standard Tomcat instance.


texts: I'm trying to get R into QGIS on my mac. I've installed the plugin "Processing R Provider plugin" but I can't "get R scripts from online scripts collection". It has something to do with the all caps "PATH" variable or something. The other similar questions like mine asked on here haven't been any helpThe first picture is what my processing toolbox looks like currently. The next two pictures are what I'd like it to look like


texts:
I am trying to determine if anyone has attempted/done directed sampling in any field. I am having a hard time trying to determine the best way to do this, but I'll try to explain.

Given a yield map like above, I want to determine areas that need sampling using statistics. Basically, I want to break down this field into areas of a certain size, most likely 75-100 ft due to memory and then determine areas that need to be sampling based on the most extreme anomalies (either very high or low) and start heading toward 0 deviation.

I need to be able to control the number of points to look at, for cost reasons.
An example is below.

I realize this has many parts, but I am only asking for direction, or examples if someone has done something similar before.


texts:
im trying to write a script that allows me to iterate through my map layers, and for each layer, iterate over each feature and add a field that gives the location of the feature using my index layer as reference. In other words, I want to know where every single feature is in relation to the index layer. The index layer is comprised of many polygons called "zones". So every feature should have a field that tells me which zone it is in. I need some help with the expression
this is what I have so far :

            import arcpy
arcpy.env.overwriteOutput=True
mxd = arcpy.mapping.MapDocument(r"C:\Users\Daimon Nurse\Desktop\Zonetrial.mxd")
df = arcpy.mapping.ListDataFrames(mxd,"Main Map")[0]
lyrList = arcpy.mapping.ListLayers(mxd,"",df)
fieldname = "Zone"
type = "TEXT"
inlayer = "layeritem2"
inlayer1 = "Zoneitem"

#for layeritem in lyrList:
    #if layeritem.name == "Zone of Interest":
        #Zones = set([r[0] for r in arcpy.da.SearchCursor(layeritem, "ZoneName")])
       #print Zones
for layeritem in lyrList:
    if layeritem.name == "Zone of Interest":
        with arcpy.da.SearchCursor(layeritem, "ZoneName") as cursor:
            for Zone in cursor:
                print Zone
                arcpy.MakeFeatureLayer_management(Zone,inlayer1)
                print Zone
                for layeritem in lyrList:
                    if layeritem.name == "AllZones_Topology":
                        print "Nothing will happen to " + layeritem.name
                    elif layeritem.name == "UWI Campus Boundary":
                        print "Nothing will happen to " + layeritem.name
                    elif layeritem.name == "Zone of Interest":
                       print "Nothing will happen to " + layeritem.name
                    elif layeritem.name == "polytrial selection 5":
                        print "Nothing will happen to " + layeritem.name
                    else:
                        print layeritem
                        print Zone
                        arcpy.AddField_management (layeritem,fieldname,type)
                        arcpy.MakeFeatureLayer_management(layeritem,inlayer)
                        #arcpy.MakeFeatureLayer_management(Zone,inlayer1)
                        arcpy.SelectLayerByLocation_management(inlayer,"intersect",inlayer1)
                        arcpy.CalculateField_management(inlayer,fieldname,Zone)

print "Script Finished!"


texts: Styles from Geomedia (mdb with style library) to be used in qgis?

texts:
We are trying to move some projects from Geomedia Professional to QGIS. Since we have large style libraries (stored as mdb) we are looking for a way to convert the styles to a format that QGIS can read. There is no export functionality in geomedia. Export in Access from mdb to xml failed when importing in stylemanager in QGIS. Can anyone help?

(Windows 7; Geomedia Professional 13; QGIS 1.8.0)


texts:
I have some problems with PyQGIS while generating a layout.

Here is a minimal (non-)working example in which I first load some layers with their styles. Then I would like to export them using a layout exporter but I only have blank files.

EDIT: simpler code

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import logging
import os

from PyQt5.QtCore import QSize
from qgis.core import *

logger = logging.getLogger(__name__)

# supply path to qgis install location
QgsApplication.setPrefixPath('/usr', True)

# create a reference to the QgsApplication, setting the
# second argument to False disables the GUI
qgs = QgsApplication([], False)

# load providers
qgs.initQgis()

# Write your code here to load some layers, use processing
# algorithms, etc.

project = QgsProject.instance()  # type:QgsProject
project.setCrs(QgsCoordinateReferenceSystem('EPSG:4326'))

dbpath = os.path.abspath(os.path.join('MWE.sqlite'))
dirpath = os.path.dirname(dbpath)

# Load cities layer
cities_style = os.path.join(dirpath, 'Style_Cities.qml')
display_name = 'Communes'
table = 'cities'
uri = dbpath + '|layername={}'.format(table)
cities = QgsVectorLayer(uri, display_name, 'ogr')
if not cities.isValid():
    msg = 'The loading of the layer named {!r} from the uri: {!r} failed.'.format(table, uri)
    logger.critical(msg)
    raise RuntimeError(msg)

# Add the layer without showing it
project.addMapLayer(cities)  # ,False)

# Style the layer
cities.loadNamedStyle(cities_style)

# Load additional layers
# ...
#

#
# Layout
#
layout = QgsPrintLayout(project)
layout.initializeDefaults()
layout.setName('console')
layout_page = QgsLayoutItemPage(layout)
layout_page.setPageSize(QgsLayoutSize(210, 297, QgsUnitTypes.LayoutMillimeters))
layout_page.decodePageOrientation('Portrait')

extent = cities.extent()  # type:QgsRectangle
extent.grow(0.1)

layout_map = QgsLayoutItemMap(layout)
layers = project.mapLayers()
layers_list = layers.values()
layout_map.setLayers(layers_list)
crs = QgsCoordinateReferenceSystem('EPSG:3857')
layout_map.setCrs(crs)

layout_map.setExtent(extent)
layout_map.refresh()
# layout.addLayoutItem(layout_map) # Cause the error on PNG export
layout.setReferenceMap(layout_map)

layout_exporter = QgsLayoutExporter(layout)
image_export_settings = QgsLayoutExporter.ImageExportSettings()
image_export_settings.dpi = 200
image_export_settings.imageSize = QSize(2000, 2000)  # pixels

res = layout_exporter.exportToImage(os.path.abspath('toto.png'), image_export_settings)
if res != QgsLayoutExporter.Success:
    raise RuntimeError()

pdf_export_settings = QgsLayoutExporter.PdfExportSettings()
pdf_export_settings.dpi = 200
res = layout_exporter.exportToPdf(os.path.abspath("toto.pdf"),     pdf_export_settings)
if res != QgsLayoutExporter.Success:
    raise RuntimeError()

# Exit
qgs.exitQgis()

I have the following error for the PNG export:

ERROR 6: The PNG driver does not support update access to existing datasets.

EDIT: if I comment the line layout.addLayoutItem(layout_map), the error disappears. It makes sense because I already provide layout to the constructor of the layout_map object. The link is already existing.

Nevertheless, I still get blank files.

EDIT: here is a link to download a valid database and style files to make the MWE work: https://we.tl/t-pwuEm8WXCd


texts: Geonetwork - Editing metadata - ERREUR 400 : Bad Request

texts:
I'm using Geonetwork 3.8.2, and I can't edit a metadata.

Th request : http://localhost:8080/geonetwork/srv/api/records/3621/editor?&currTab=default&starteditingsession=yes&_random=9472

return an error 400 :

  <apiError>
    runtime_exception
    <description>Failed to compile stylesheet. 12 errors detected.</description>
    <message>TransformerConfigurationException</message>
  </apiError>

My config is :
Tomcat 7.0.76
Open JDK 1.8.0_232
Geonetwork 3.8.2

Has anyone ever encountered this problem?


texts: Georeferencing Vector Layer (shp) from Mollweide Projection to Custom CRS in QGIS

texts:
I'm working on a fantasy world; a map, to be more specific. I basically just started to paint a world map in Photoshop, but over time it became very big an unfortunately Photoshop is not really good with detailed world maps.

So I went out to look for another software to manage coastlines, rivers, cites, town, borders, streets and much much more. I decided to use QGIS, imported my .psd as .tif and started to vectorize it manually (yes, it took me quite long and I haven't finished yet).

The world map is a projection of a real planet. I never set up any CRS. The problem is, I never thought of any projection-issues before.
Later I decided, that the painted map is a "Mollweide" projection of my planet. (I know this will cause some distortion, but that's okay)

To make it easier, I decided the planet to be a perfect sphere with r=10567km.

Since I am a newbie on QGIS, I have three issues now:

I have no good idea, how to define a custom CRS and there is no tutorial for doing that.
I have no idea how to georeference my vector layer. I worked on a mathematical way to transform the mollweide-x,y-coordinates in QGIS to real latitude-longitude-coordinates. But I don't know how to perform the algorithm on my shapefile (I thought of exporting it to GeoJSON) and get QGIS that to display the geocoordinates using my custom CRS after transforming.
Is there a way to keep my background raster layer after transforming the vector layer so I can go on vectorizing?


texts: Esri ArcGIS shapefile to CAD - CAD file returns empty?

texts:
I'm using ArcGIS 10.6.1, and I need to transform a points shapefile into a CAD file for AutoCAD.

I used the "to CAD" transformation tool in ArcMap and everything seemed to go smoothly, but when I tried to open the CAD file in TrueView, nothing came up. It wasn't a visibility issue because the layer was visible in TrueView.

What is the problem? I should also caveat that it's my first time doing this transformation. I have plenty of experience doing conversions the other way round though.

I used a coordinate system that I've created; this custom coordinate system is what I use on the CAD files the same architects give me to line up my data with their true geographic location instead of CAD's local coordinate system

The point file's attributes are: (1) latitude (2) longitude (3) name - nothing else! How does a shapefile's attributes affect how it gets transformed into a CAD file?


texts: tablejoin in qgis 2.x adds name of table to columnheader attributetable

texts:
When joining a table to a shapefile through 'layer properties > joins' the columnheader of the attributetable shows the name of the file you joined.

eg.:
the name of the table is: STAVAZA3

the table looks like:

after joining with a shapefile, the attributetable looks like:

This didn't happen in qgis 1.8.x.

Predefined legends wouldn't work anymore because of the changed names of the columns.

Is this the standard way how qgis 2.x joins tables?  Or is it due to a property that i have to set anywhere?


texts: How to add new projection using ArcObjects of ArcGIS for Desktop?

texts: I want to use "Mars MOLA" { 86 3396000.0 3396000.0 "Mars MOLA" 0 0 0 0 0 0 0 } with my shapefiles.


texts:
If I move the layer's labels with the move label tool my label is overlapping. The labels are categorized by rule

My rule are based on two field "column1" and "column2"

I can't use only one field because I need a different color for every single rule
How I can I fix it?


texts:
im able to create a data store in geoserver but when i try to access it in the geoserver admin console i get this error :

can't load chosen data, verify that the loader module is installed (sorry for this horrible traduction because the error message is in french)

and this is the code :

   public function createGeoTiffDataStore($datastoreName, $workspaceName, $location) {
    $data = '<dataStore>
                  <name>'.htmlentities($datastoreName, ENT_COMPAT).'</name>
                  <type>GeoTIFF</type>
                  <enabled>true</enabled>
                  <connectionParameters>
                    <entry key="url">file:'.htmlentities($location, ENT_COMPAT).'</entry>
                    <entry key="namespace">'.htmlentities($workspaceName, ENT_COMPAT).'</entry>
                  </connectionParameters>
            </dataStore>';
    return $this->runApi('workspaces/' . urlencode($workspaceName). '/datastores', 'POST', $data);
}

private function runApi($apiPath, $method = 'GET', $data = '', $contentType = 'text/xml') {
    $ch = curl_init();
    curl_setopt($ch, CURLOPT_URL, $this->serverUrl.'rest/'.$apiPath);
    curl_setopt($ch, CURLOPT_USERPWD, $this->username.":".$this->password);
    if ($method == 'POST') {
        curl_setopt($ch, CURLOPT_POST, true);
        curl_setopt($ch, CURLOPT_POSTFIELDS, $data);
    } else if ($method == 'DELETE' || $method == 'PUT') {
        curl_setopt($ch, CURLOPT_CUSTOMREQUEST, $method);
    }

    if ($data != '') {
        curl_setopt($ch, CURLOPT_HTTPHEADER,
            array("Content-Type: $contentType",
                'Content-Length: '.strlen($data))
        );
    }

    curl_setopt($ch, CURLOPT_HEADER, false);
    curl_setopt($ch, CURLOPT_RETURNTRANSFER, true);
    $rslt = curl_exec($ch);
    $info = curl_getinfo($ch);

    if ($info['http_code'] == 401) {
        return 'Access denied. Check login credentials.';
    } else {
        return $rslt;
    }
}


texts: Persistence in ArcObjects 10.4 back to IPersistVariant?

texts:
In the 9.x versions of ArcGIS Desktop the way to persist data in the MXD was through implementing the IPersistVariant interface. In version 10.0 the advised way to persist data was use of the PersistenceHelperClass.

Now in 10.5 the IPersistVariant is advised again, there is no mention of the PersistenceHelperClass anymore and my code fails. Is this really happening?


texts: Save WFS layers multiple records with geoext editor grid

texts:
How can I save multiple records after modification with geoext editor grid?
The saving is working fine if I modify only one of the records, but it fails when I try to save multiple modifications.

I used this geoext solution to modify the attributes of my wfs layer:
How to edit feature attributes with openlayers?.


texts: QGIS Server: list metadata of all available services and layers

texts:
I would like to autogenerate a landing page for a QGIS Server where all available services and layers with the corresponding metadata is available.

It is meant as documentation for within an organisation. I am not looking for a webclient. I am looking for something similar Esri provides for ArcGIS servers.

Is this possible?


texts: Filtering Sentinel2 image collection to get least cloudy image for entire ROI

texts:
I am trying to obtain the least cloudy image (or images) across my entire ROI. Though I do not get any errors, the true color composite is not being mapped for me to check if I have indeed sorted and used the least cloudy image in the entire Sentinel 2 image collection in my ROI and within a date range.

/////////////////////////////////////1. Parameters
//Load asset
var asset= ee.FeatureCollection('users/myfolder/asset'); //you can use any asset you have, eg consider a country like US. I cannot make my asset public
Map.addLayer(asset);

//Load Sentinel 2 image collection
var sentinel2 = ee.ImageCollection('COPERNICUS/S2');

//Function for masking clouds
// http://xzsunbest.tk/2018/07/05/MaskingOutCloudInSentinel2WithGEE/

// Bits 10 and 11 are clouds and cirrus, respectively.
var cloudBitMask = ee.Number(2).pow(10).int();
print (cloudBitMask); //1024
var cirrusBitMask = ee.Number(2).pow(11).int();
print (cirrusBitMask); //2048

var cloudmask = function (image) {
var qa = image.select('QA60');
// Both flags should be set to zero, indicating clear conditions.
var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
         qa.bitwiseAnd(cirrusBitMask).eq(0));
return image.updateMask(mask);
 }

// Create clipping function
var clipper = function(image){
return image.clip(asset);
  };

/////////////////////////////////////2. Application
var clipped = sentinel2.map(clipper);
//var filtered = clipped.filterDate('2015-06-01','2018-08-30')
//                  .filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE',2))
//        filtering to images with cloud cover of 2%. I do not want to
//filter.
//print (clipped); //Image Collection- will not print because there are >5000
//elements
var count = clipped.size();
print('Count: ', count); //8168695 across Sentinel2 in ROI

var filtered_date = clipped.filterDate('2015-01-01','2015-12-31');
var count_filtered = filtered_date.size();
print('Filtered (by date) Image Count: ', count_filtered);//259629 reduced
because of date range

var images = ee.Image(filtered_date.sort('CLOUDY_PIXEL_PERCENTAGE').first())
print('Least cloudy image: ', images);
Map.addLayer(images, {bands:['B4','B3','B2']}, 'true color composite');


texts:
I'm trying to change OpenLayers.Lang file to Turkish. I succeeded to some extent except for encoding. when I use şiüçöİı chars that are particular to Turkish, they are replaced with some strange characters.

Is there a way to handle this problem?


texts: Finding the date when a specific polygon was created using ArcGIS Desktop

texts:
I have a shapefile of polygons. Polygons were created in different times (digitized on raster image). I would like to determine the date that specific polygon was created.
Is there any tool that can help to find this out in ArcGIS Desktop?


texts:
I need done seems to be a simple task, but I cannot find out how.

I want to create grid lines for a specific area at a specific angle (which is the easy part).

Now I want on every gridline to be shown the lowest elevation and highest, and by multitude of x every point on the line where the elevation reached the higher mark.

I'm a newbie and self taught using google earth and global mapper.


texts: Grouping layers by attribute value and finding nearest neighbour using QGIS?

texts:
I am new to QGIS.

I have two layers (layer_a and layer_b) consisting of two sets of points.

If I want to find the nearest point in layer_b for each point in layer_a, I can just use the Distance matrix calculator that comes with the graphical interface in QGIS.

My issue is slightly more complicated. I want to group both layers by an attribute value, before I calculate the distance between points.

Translation: if layer_a has a bunch of points to libraries around the world and layer_b has a bunch of points to big cities around the world, and both layers have an attribute where the value tells me which country each point is in, then I want to find the shortest distance between a library point in a specific country and a big city point in the same country.

I could accomplish this by splitting the layers by country and then repeating the analysis on each layer (i.e. distance between library points in Germany and big city points in Germany), using the graphical interface in QGIS, but this is tedious. Also, in reality I will have thousands of groups.

How do I do this in QGIS?


texts: Write Nominatim-like geocode query using only PostGIS functions?

texts:
I'm new to PostGIS and was anticipating being able to write geocode queries using only the functions provided by PostGIS.  However, once I loaded my database using osm2pgsql and discovered tables like "point," "polygon," etc., it wasn't at all apparent how I could create the equivalent of a Nominatim query such as this one for a street address in Cologne, Germany:

# This returns "lat":"50.9412187","lon":"6.95511930515345"
http://nominatim.openstreetmap.org/search?q=margarethenkloster+5,+50667,+cologne,+de

I know that latitude/longitude data is kept in the nodes table but the columns found in the point, line, and polygon tables don't seem to lend themselves to this type of query.

Is it possible to write a SQL query against a PostgreSQL/PostGIS database that will mimic the above call to the Nominatim web service or must I use Nominatim?

Thanks.


texts: How to find dead-ends in road network using ArcGIS for Desktop or GRASS?

texts:
I'm trying to remove all the dead-end roads in a road dataset (lines) that are less than 100m long.

How do I identify roads that are dead-ends (no exit)?

The 100m part is easy.

I know how to use ArcGIS for Desktop proficiently and can use GRASS.


texts: Clip a raster by vector in QGIS while keep raster cells that do not touch the vector boundary

texts:
I am using QGIS 2.18.11 with the Semi - Automatic Classification Plugin (SCP). I use Sentinel 2a as well as Landsat 8 Scenes and run some Classifications with SCP.

Inside my Sentinel 2a scene I have several areas I want to analyse. Therefore I want to clip the (classification)raster by a vector (I created upfront). This works perfect with "clip multiple rasters" by using the vector ... BUT: I need just the pixels inside the vector boundaries which do not have contact with this vector boundary.

I tried an option but it does not succeed:

I transformed the clipped classification raster to vector and clipped the two vector data (boundary and the tranformed classifiaction raster). That worked.

Then I tried to transform this new vector to a raster and hoped that there will be no pixels lying on the boundary but I just get my old classification raster back.

Any comments or suggestions?

Pixel goes beyond the boundary

This looks like as it should be


texts: Configuring URL link for popup window directly in qgis2web plugin?

texts:
I am using the qgis2web plugin (the leaflet mode) and everything works fine. I created few objects in qgis as a points in shapefile layer on the map.

What i need now is to add some url link to each point's popup window (now i can only put there id number of the point, but not the url link).
I can do it afterwards by changing html code of the final website, which is generated by the qgis2web plugin, but is there any way to do it right in qgis? It would be much more easier to do it right in qgis.

If is it possible, can i add different url link to each point of one layer or the url link have to be always the same for each object of the layer? What i finally need is a different url link in the popup window of each point, which opens new small window with another website.


texts: What is the complete list of values of ST_GeometryType() in PostGIS?

texts:
What is the complete list of return values for PostGIS ST_GeometryType()?

The documentation linked above only says:

  Returns the type of the geometry as a string. EG: 'ST_Linestring',
  'ST_Polygon','ST_MultiPolygon' etc.

I'm trying to determine the type of a geometry in a function, but couldn't find the list of all cases. (This is with PostGIS 2.5/2.4)


texts: Link Excel Table to Access/Arcmap Personal Geodatabase Table

texts:
My overall goal:

My coworkers have an excel table which they update/use frequently.  The table is not in a format that arcmap can decipher.  Therefore, I've created a locked sheet (in the format I need) within the excel table, which points to the data my coworkers use.  Periodically throughout the year I need to create a map with features that reflect the changes made in said excel table. My thinking is to store my own table in a geodatabase along with my feature class. I want to be able to easily update my table via a linked excel sheet within said gdb. So ultimately for every new map, I would like to just do a quick query in Access to update my table instead of manually doing it every time.

What I've done so far:

I know arcmap does not recognize the linked table alone, so I must add the linked data to a table, hence this next step:
I created a personal geodatabase (.mdb) in arcmap and created a new table within the gdb which includes all the fields I need.  When I open my gdb in Access 2013, I add my linked excel sheet.  (I want to keep the linked table so it automatically updates when the excel table is updated).

After linking my excel sheet, I want to the append the values of the linked excel sheet to the table I created in my gdb.

I've tried using the Append Query, but nothing happens.  The linked excel sheet values do not append to the table I created in my gdb.  I also don't see an error messages.

Is there another way to do this?  Am I doing something wrong?


texts:
I have inverted the color ramp to display land surface elevation, now when I add a legend in print composer it shows the lowest value at the top and the highest at the bottom. See below:

It is more logical to invert that too (so +10 is at the top and -10 is at the bottom). Is this possible?


texts: Getting road name, road type and max speed for many GPS coordinates

texts:
I have compiled a local instance of the overpass api, because I need to make a lot of "around" queries (over 6 million). The solution I'm using right now is to call bin/osm3s_query for every single query. This is very slow, because osm3s_query has to start, answer the query and shut down for every single query.

The solution in Query one overpass api process multiple times (which question is very similar to mine but not the same) suggest just grabbing the bounding box and then searching the return value of that, but for "around" queries that doesn't help because I have no perfect matches for the GPS coordinates that gets me. Is there a way to submit multiple queries to the overpass API? Or should I be using something else entirely?

edit: it currently takes a little under half a second to get a single point. Even if I cut down my dataset an order of magnitude that's still 3.5 days non-stop.


texts: Filtering date from images of Earth Engine JRC database

texts:
I am new in Google Earth Engine. I used JRC image collection which is an image of global surface water from 1984-2019.

var gsw = ee.Image('JRC/GSW1_1/MonthlyHistory');

The problem is I cannot filter date from the image.
I want to filter the date from 1994 to 2003.

Is it possible to filterDate from the image or filterDaate is only applicable for images collection?


texts:
QGIS 3.0 is released. However, in the download page, there does not seem to be an entry (for Ubuntu (16.04) at least):

  Latest Release    2.18.x Las Palmas (new LTR)

  https://qgis.org/debian

  ...

  Long Term Release Repository  2.14.x Essen (previous LTR)

  ...

  Development Version   2.99.x master [4]

  https://qgis.org/debian-nightly

  ...

Which repository should be used to get the new QGIS 3.0?


texts: How to create a QGIS road style for OSM data?

texts:
I am trying to use qgis to visualize the shortest path algorithm's result. I managed to import the osm map into qgis

But I think it looks kind of different. I was referring to this link and she has this nice map like this

So I tried to import the database creating using osm2po. And this is what I got

This is still different from the one seen at link. I know the region is different, but still

I followed the same step as defined in link link. I also used the same file for austria and imported it into my database. I changed the style as well and this is what I got.

Now the roads and highways in the link link look nice and fine. However, in my case they look kind of overlapped. Why is it so?


texts: How can I set the project CRS to an existing User Defined Coordinate System with python in QGIS?

texts:
I want to set the project CRS with python in QGIS. If I set it to an EPSG code, it works as expected:

canvas = iface.mapCanvas()
selectedcrs="EPSG:4326"
target_crs = QgsCoordinateReferenceSystem()
target_crs.createFromUserInput(selectedcrs)
canvas.setDestinationCrs(target_crs)

But when I want to set it to a User Defined Coordinate System, the project CRS is made empty, because it's code is not in the predefined CRS list in QGIS, but it's in the User Defined Coordinate Systems list in QGIS.

canvas = iface.mapCanvas()
selectedcrs="USER:100001"
target_crs = QgsCoordinateReferenceSystem()
target_crs.createFromUserInput(selectedcrs)
canvas.setDestinationCrs(target_crs)

Is there a way to set the project CRS to an existing User Defined Coordinate System? Or is there a way to get the definition of that existing User Defined Coordinate System with Python?

EDIT:
To be clear: I don't want to change the CRS of a layer. I want to change the CRS of the project.

EDIT:
I can select "USER:100001" and set it from the Project Properties, but I want to do that with Python.

EDIT:
In my full script the "USER:100001" comes from the Coordinate Reference System Selector (which also lists the User Defined Coordinate Systems) and  use it with this code:

projSelector = QgsGenericProjectionSelector()
projSelector.exec_()
projSelector.selectedCrsId()
selectedcrs=projSelector.selectedAuthId()

The selectedcrs variable is stored as a setting and later used by the script I originally posted above.


texts: Overwriting Hosted Feature Service via REST

texts:
A quick skills background... Not a full time programmer, proficient enough in Python to get by but VERY new to REST API scripting. So forgive me if there is an easy answer, I am missing something simple or a different approach is warranted.

I am trying to update an existing hosted feature service by using the 'Publish Item' operation using the overwrite=true parameter. (API Reference: ArcGIS REST API - Publish Item).

The following function (minus overwrite param) works perfectly to publish the feature service from an existing geojson item I have already added to Portal.  When I add the overwrite parameter, I get the message that a service with that name already exists. Seems like it checks the service name before recognizing the overwrite parameter.  Service is not delete protected.  Seems simple and straightforward adding of a parameter to the REST call.  Resulting URLs from the original publishing POST call and the update call are the same except for '&overwrite=true' is added to the URL.

What am I missing?

Response Message:
"Publish error for item 'CMT_Test3.geojson'. Service name 'CMT_Test3' already exists"

def PublishService(Short, itemID, username, token, itemName, publishParams):
    publishURL =   'http://{}/portal/sharing/rest/content/users/{}/publish'.format(Short, username)
    query_dict = {
        'itemID': itemID,
        'filetype': 'GeoJson',
        'f': 'json',
        'token': token,
        'overwrite': 'true',
        'publishParameters':publishParams}

    query_dict['publishParameters'] = json.dumps(query_dict['publishParameters'], sort_keys=False)
    request = requests.post(publishURL, data=query_dict)
    return request.json()


texts:
I open a shapefile with the rgdal lib:

foo <- readOGR(dsn="myfolder",layer="myshapefileWithoutEnding")
foo.df <- as(foo, "data.frame")

Then I add a column to the data frame with

foo.df$MY_VALUE <- myDataFrameWithOneColumn

The I check if it looks fine with print(foo.df) and it does, foo.df shows my values as a new column

Then I write the shapefile

writeOGR(foo, ".", myfilename, driver="ESRI Shapefile")

It saves 4 files, everything looks fine.

But when I open the shapefile in a viewer, I can't see my added column. Everything looks like before.

Any idea why?


texts: How to identify proper CRS for new point layer?

texts:
I am having trouble aligning a new point layer. I created a new layer and uploaded a CSV UTF-8 file with the projection in my project but can't seem to get it to align properly. I am working with
EPSG:6499 - NAD83(2011) / Michigan South (ft) - Projected
Here are a few points I cant seem to align.
What am I missing here?


texts: Double transformation of marker in Carto CSS

texts:
I need to show on the map the direction of the road traffic with arrows. For this purpose I use markers for Carto CSS in TileMill.

If the road has in OSM tag oneway="yes" or oneway="1", all is ok, but sometimes the tag is oneway="-1" that means my direction arrows should be rotated to 180 degrees. It would not be a big problem (marker-transform: "rotate(180,0,0)";) but at the same time I also need to scale my markers in different zooms (for ex., marker-transform: "scale(0.5)";).

I have no idea how to apply both transformations at the same time - I change their positions in the code but it shows me rather only first or second transformation.


texts:
I need some help writing the sql for returning the "MEAN" value from a raster(in postgis db) that intersects a defined(inquery) polygon(in postgis db).

The polygon layer is called: public.farmpaddocks
The raster layer is called:  public.l8_9385_20131118ndvi_vg94
Both are SRID :3111

I have a field in the public.farmpaddocks table that I would like to filter the result on.
Field: farmpaddocks.farmid = '511'.

I assume i need to use the ST_Summarystats function? along with intersection and clip ??
This is my first attempt at an intersect. I get no error or result from the pgadmin output pane?? Not sure whats going on.

sql:

SELECT
l8_9385_20131118ndvi_vg94.rast,
farmpaddocks.farmid
FROM
public.l8_9385_20131118ndvi_vg94,
public.farmpaddocks
WHERE
farmpaddocks.farmid='511' AND ST_Intersects(l8_9385_20131118ndvi_vg94.rast,
farmpaddocks.the_geom);


texts: Create a geography/geometry column from x and y fields - SQL Server 2008

texts:
I'm new in SQL Server 2008 and I hope you will understand my question/need.

Thus, I have a table which contains 3 fields (Name, Lat and Long) in my data base (spatial). I want to create a geometry/geography column based on those fields (Lat and Long) but unfortunately without any success.

My question is: How can I manage to do that?


texts: How to reverse quantile color palette with CARTO VL

texts:
Using CARTO VL (which is awesome by the way), I have made the following groundwater harvest potential map that is styled using 10 quantiles.

Although the "temps" palette works well, is it possible to reverse the colors? I read a thread on GitHub that mentioned using "inverse" like this:

color: ramp(globalQuantiles($sa_m3_km2, 10), inverse(temps))

But this have never worked for me. Is inverting a quantile palette possible?


texts: Polygon not rendering on my map OpenLayers3

texts:
EDIT: I added a fiddle here: https://jsfiddle.net/gnsmq2t9/

I am new here, and new to GIS I am building an app, and I want to put a polygon around different sections that will allow the user to click and zoom into that section. I am having problems getting my polygon to appear, I am not getting any console errors.

Here is my source/feature

    var myOtherVectorSource = new ol.source.Vector({});
    var polyCoords = [
        ol.proj.transform([-121.901559, 37.417891], 'EPSG:4326', 'EPSG:3857'),
        ol.proj.transform([-121.901139, 37.416427], 'EPSG:4326', 'EPSG:3857'),
        ol.proj.transform([-121.899634, 37.416678], 'EPSG:4326', 'EPSG:3857'),
        ol.proj.transform([-121.900312, 37.418398], 'EPSG:4326', 'EPSG:3857'),
        ol.proj.transform([-121.901559, 37.417891], 'EPSG:4326', 'EPSG:3857')
    ];
    var pLotSection1 = new ol.geom.Polygon([
        [
            polyCoords
        ]
    ]);
    var pLotSection1Feature = new ol.Feature({
        name: 'Parking Lot One',
        geometry: pLotSection1
    });
    myOtherVectorSource.addFeature(pLotSection1Feature);

Here is my layer

var myOtherVectorLayer = new ol.layer.Vector({
        source: myOtherVectorSource
    });

And my map

myLayers = [myStaticImageLayer, myVectorLayer, myOtherVectorLayer];

    //set the view, center, and zoom level
    var myView = new ol.View({
        center: ol.proj.transform([-121.896436, 37.413208], 'EPSG:4326', 'EPSG:3857'),
        zoom: 19,
        enableRotation: false,
        constrainRotation: 0
    });
    var map = new ol.Map({
        target: 'mymap',
        layers: myLayers,
        view: myView
    });

I have another vector/feature layer that I am plotting similarly, and that is showing up perfect. Please let me know if I posted incorrectly. Like I said, am new here and to GIS.


texts: Filter option in Semi-Automatic classification Plugin, QGIS

texts: I usually download Sentinel-2 scenes using SCP plugin in QGIS, but I don't know how to use the Filter option in this process.


texts: When editing layer, cannot save edits for lines or polygons

texts:
I am trying to create a shapefile to draw a line.I can toggle edit, but the line does not correspond to the shapefile (when I unstick the layer, the line remains) and the save option is blacked out. The line is not 'selectable' either - once Ihave drawn it I cannot remove it.

I am using QGIS version 2.18.3. Anyone have any idea why?


texts:
I've been trying to draw 3D flight path in Cesium without any success. I followed the answer from Drawing flying/moving line on Cesium?

However, this only renders a segment of flight path floating in the air. What I am trying to do is to render the flight path segment that can be projected/extruded to the ground as in the image below. I looked at the Cesium API document for PathGraphic options https://cesiumjs.org/refdoc/. This API only has leadTime and trailTime.

Does anyone know how Cesium can render a flight path like this?


texts: How to change the attribute values of the feature using the custom feature form

texts:
When i am using the following command it opens the feature form with all the fields:

qgis.utils.iface.openFeatureForm(qgis.utils.iface.activeLayer(), f)

But i want to open a feature form with the specified field names. How can i do that using python programming.


texts: How to determine elevation drop along stream with ArcGIS Desktop?

texts:
I need to calculate(estimate) hydraulic head(elevation drop) along a single river stream in order to find potential site for small hydropower plant.

I have a DEM with 20m resolution and raster and vector format of the stream.

Can anybody help me how can I do that with ArcGIS tools or in some other way?


texts: How to deploy ArcGIS Engine (Standalone) with License Agreement?

texts:
I've completed a standalone app in ArcGIS Engine and in the process of deploying it, but I have a license agreement and supposed to add a licenseagreement.rtf file format in the properties menu while setting up in visual studio 2010.

Am I to create this file myself or am I able to find it?


texts: QGIS: Placing points on a shapefile in relation to their distance to other points on that shapefile

texts:
I'm attempting to add UK railway signal-boxes to QGIS to extract their coordinates - there's over 1000 of them.
I have:

The entire UK railway network shapefile
The coordinates of train stations on the railways
The distances of the signal-boxes from the closest train station

I was planning to obtain the signal-box locations based on their distances from other train stations.
However, I don't know a direct way of splitting the railway network into smaller segments in order to use methods like Linear Referencing. I was also hoping, if it's possible to segment the network, to name each line segment and add the signal-boxes based on this name.
Any suggestions?


texts: Can the SAGA algorithm "Local Minima and Maxima" not find any minima and maxima if the spatial resolution is too high?

texts: I work with a 5cm high resolution CHM to determine tree tops. However, the algorithm does not work. If I resample my data set to 25cm, I also get a result.


texts: Polygon geometry suddenly disappeared/deleted from shapefile in QGIS

texts:
I was editing a polygon shapefile in QGIS (3.0.1) today when suddenly 966 out of 5382 polygons disappeared from the map. I had been changing attribute information in the table, cutting polygons, etc as I was going along and saving my shapefile regularly. The polygons are still in the attribute table (if I select all my polygons using a spatial selection and then inverse the selection, the lines with the missing polygons show up) but an error message (zoom failed, feature has no geometry, "le zoom vers l'entité a échouché, l'entité n'a pas de géométrie") appears when I attempt to find those polygons. There was no spatial clumping to the lost polygons nor a common attribute.

Comparison of yesterday's file (purple underlay) with today's corrupted file (black hatched)

I can't share the file because it is proprietary so I realize that it is limited what can be troubleshooted just on my description. I have the back-up in our server from end of last workday so I will have only lost one day of work but it is still frustrating.


texts: Show groups of coordinates on a map, to visually find coordinates that belong to the wrong group

texts:
I have a (very) long list of casinos. Each casino has:

latitude
longitude
boss (a single boss runs from 1 to 50 casinos)

Sample:

40.642 -73.997 Freddy Queens
40.680 -73.985 Freddy Queens
40.697 -73.949 Freddy Queens
40.651 -73.968 Freddy Queens
40.772 -73.981 Manhattan Michael
40.813 -73.939 Manhattan Michael
40.755 -74.000 Manhattan Michael
40.781 -73.964 Manhattan Michael
40.721 -74.062 Nicole Je. Polizzi
40.750 -74.042 Nicole Je. Polizzi
40.694 -74.090 Nicole Je. Polizzi
[... 70,000 casinos for 27,000 bosses ...]

From this data I want to generate a map like this:

Ideally one random color per boss. In my world the casinos of each boss are supposed to not "overlap" too much with other bosses, so my final goal is to use this map to find anomalies (casinos that should be transferred to a different boss because they are clearly in that boss' "territory"). Boss "territories" do not match any official geographical subdivisions, for all that matters the same problem could be set on Mars (Earth map background would be nice though).

Is there a webapp to which I could feed this data (CSV) and that would show this kind of map? It can look different, for instance monochrome areas would be OK too, as long as it can be used to detect anomalies. A webapp would be the best as many GIS-illiterate people need to quickly check it every once in a while.

I have tried geojson.io but it can not colour points, all points appear equal.

I have tried Mapbox, but it does not seem to support metadata (such as boss) import, and even less modifying a marker's style based on that metadata.


texts: Suggest an algorithm to measure how a node is connected to road network

texts:
I try to measure how well a node is connected to the rest of the network. If a city has a distant suburb connected with one sole road, it's connectedness is, say 1. If it's connected by 2 roads, it's 2, and so on. A point in the middle of the city is connected infinitely.

What are plugins or algorithms available?

I personally prefer open source solutions (QGIS, PostGIS, or python-based), but other options are also welcome.


texts:
In the arcpy window in arcgis, I'm trying to:

iterate through each layer in an mxd
clip the layer to a shapefile (also in the mxd)
export each clipped layer to a geodatabase, with the name of each new clipped feature class being the same as the input layer.

My code:

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, '')[0]
import os
workSpace = r"J:\foo\bar\test.gdb"
for lyr in arcpy.mapping.ListLayers(mxd, '', df):
    fullPath = os.path.join(workSpace, lyr)
    arcpy.Clip_analysis(lyr, "clipping_shape", fullPath)

Error message:

TypeError: object of type 'Layer' has no len()

Edited to add: above error message appears even without the clip_analysis portion of the code, so presumably something is off with the fullPath variable.

I've read this answer but I don't know how to use this portion of the code when populating a geodatabase instead of using shapefiles.

"{}.shp".format(lyr)


texts: Converting historical (1992/2000/2006) TIGER/Line files to modern shapefile formats

texts:
I'm trying to get some historical geographic information from Wisconsin and have been pointed the way of the TIGER/Line Shapefiles. Unfortunately, they're in a strange format (see, e.g., the files for Adams County from 1992 [.zip download link]), with file types like .F5A (a set of files of type .F5{x} where {x} is one of A, I, P, R, 1, 2, 3, 4, 5, 6, 7, 8).

This is quite foreign to me, but I've seen it suggested that it may be possible to convert these to modern shapefiles.

Is there a robust way of doing so? (It goes without saying I prefer a free option)

Alternatively, how could I open the file in a common GIS program (QGIS, SAGA, or R), whence I could convert it within those programs (or at least view the polygon boundaries, which is a step in the right direction)?


texts:
SQL error:

  ERROR:  function point(numeric) does not exist LINE 4:
  (ST_GeometryFromText(Point(53.34972 -6.26025)), 4326), 1000)
                               ^ HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

In statement:

SELECT o.name
FROM osm_point o
WHERE ST_DWITHIN(o.way,
(ST_GeometryFromText(Point(53.34972 -6.26025)), 4326), 1000)


texts:
I am new to using PostgreSQL/PostGIS.

I have two tables that I've imported from two shapefiles. One table is based on points called "shape_points" and the other based on polygons called "shape_polygons".

Now I want to get the bounding box of all records, so I use the next query for the table "shape_points":

`SELECT box2d(geom) FROM shape_points;`

AND this for the table "shape_polygons":

`SELECT box2d(geom) FROM shape_polygons;`

At this point both querys works perfect.

So the big problem is when I do a UNION of the two tables:

`SELECT box2d(geom) FROM shape_points
  UNION
 SELECT box2d(geom) FROM shape_polygons`

I get the next error:

"ERROR:  could not identify an equality operator for type box2d"

What is the problem?


texts:
Trying to login into my AGS server 10.1 via server manager throws an error
see atatched pic


texts: Export footprints of an unmanaged raster catalog with arcpy

texts:
Is there a way to access the "export footprints" function for an unmanaged raster catalog with arcpy? It can be done through ArcMap when a raster catalog is added to the map as seen below.

The context is that I am working on a python script to automate the process of adding rasters to an unmanaged catalog, and then exporting the footprints. The raster catalog is in a file geodatabase.


texts: How can I solve this error from TIFF file?

texts:
I used google colab(python) to read a tif file with Tiff file. But it shows:

TiffPage 0: <COMPRESSION.LZW: 5> requires the 'imagecodecs' package.

This is my code:
!pip install tifffile
!pip install imagecodecs
import numpy
import tifffile as tiff
img = tiff.imread('/content/gdrive/My Drive/My Driver/imageToDriveExample.tif')
matrix = img[:,:,:]
print(matrix)


texts:
My current Project is to "map" a turntable in qgis. I know, very unlikely but that's what my teacher wants.

My teacher also wants me to measure something in this project (Very specific, isn't it?). And that's where the problem starts.
Qgis wants me to define geographical coordinates or something like that. The problem is, that i don't have any coordinates because it's a picture of a turntable and not an actual map.

Any ideas?


texts: Can I customize a Google Maps layer in terms of content?

texts:
I don't want all features that come with the Google street layer. I just want basic infrastructure like roads, national parks, cities, streets with their names, etc. Everything that is public is OK. I don't want private companies and shops to appear on the map.

I am using OpenLayers to display the Google Maps street layer. So maybe it's possible to choose what to display by using the Google Maps API instead of OpenLayers..


texts: Can't compile python plugin in QGIS (Windows)

texts:
I'm just installed QGIS 2.2 and created a new plugin from Plugin Builder. When I try to compile the plugin using:

C:\Program Files\QGIS Valmiera\bin>pyuic4 -o C:\Program Files\QGIS Valmiera\bin\
myplugin\ui_myplugin.py C:\Program Files\QGIS Valmiera\bin\myplugin\ui_
myplugin.ui

I received the error message "Error: one input ui-file must be specified".

Someone could help me solve this issue? I'm running over Windows 8 OS.

Thanks,
Mateus


texts: WFS and WMS query SRS parameter and BBOX usage

texts:
WMS and WFS queries needs SRS and BBOX parameters.

http://example.com/geoserver/wfs?
  service=WFS&
  version=2.0.0&
  request=GetFeature&
  typeNames=namespace:featuretype&
  srsName=CRS&
  bbox=a1,b1,a2,b2&

For example I have GeoServer to publish WFS. My data SRS is EPSG:4326 and my map SRS is EPSG:3857. So I send a WFS query.

http://example.com/geoserver/wfs?
  service=WFS&
  version=2.0.0&
  request=GetFeature&
  typeNames=namespace:featuretype&
  srsName=EPSG:3857&
  bbox=978393.9621,2680799.4560,-430493.3433,3326539.4710&

Does GeoServer or ArcGIS server give an error? Because my data has EPSG:4326 . Or does GeoServer convert these coordinates internally and find data and send me as EPSG:3857


texts: Managing error with GPS routes (theoretical framework?)

texts:
I'm looking for the appropriate theoretical framework or speciality to help me deal with understanding how to deal with the errors that the GPS system has - especially when dealing with routes.

Fundamentally, I'm looking for the requirements on the data and any algorithms to use to be able to establish the length of a trail. The answer needs to be trustworthy.

A friend of mine was the race director of a race which was billed as 160km but the Garmin watches everybody has makes it more like 190km+. It caused quite some grief at the finish line, let me tell you!

So my friend went back to the course with various GPS devices in order to remap it and the results are interesting.

Using a handheld Garmin Oregon 300 she got 33.7km for one leg. For the same leg on a wrist watch Garmin Forerunner 310xt it came out to 38.3km.

When I got the data from the Oregon it was obvious that it was only recording data every 90 seconds or so. The Forerunner does it every couple of seconds.

When I plotted the data from the Oregon I could see that it got confused by some switchbacks and put a line straight through them and a curve was made a little less.

However, I muse that the difference in the recording frequency is much of the explanation. i.e. by recording every couple of seconds the Forerunner is closer to the real route. However, there will be an amount of error because of the way GPS works. If the points recorded are spread around the real route randomly (because of the error) then the total distance will be larger than the real route. (A wiggle line going either side of a straight line is longer than the straight line).

So, my questions:
1. Are there any techniques I can use on a single dataset to reduce the error in a valid way?
2. Does my theory about the difference in recording frequency hold water?
3. If I have multiple recordings of the same routes are there any valid techniques to combine them to get closer to the real route?

As I say, I don't really know what to search for to find any useful science about this. I'm looking for ways to establish just how long a given piece of trail is and it is very important to people. An extra 30km in a race is an extra 5+ hours we weren't expecting.

Thanks for any advice you can give.


texts: Unmet dependencies installing QGIS 2.0 on Ubuntu

texts:
I have followed the instructions-a link to install QGIS 2.0. Everything worked fine until the punch line, installing QGIS ( I do have Ubuntu 13.10);

sudo apt-get install qgis python-qgis qgis-plugin-grass
eyalak@eyalak-ThinkPad-X220:~$ sudo apt-get install qgis python-qgis  qgis-plugin-grass
Reading package lists... Done
Building dependency tree
Reading state information... Done
Some packages could not be installed. This may mean that you have
requested an impossible situation or if you are using the unstable
distribution that some required packages have not yet been created
or been moved out of Incoming.
The following information may help to resolve the situation:

The following packages have unmet dependencies:
 python-qgis : Depends: python-qgis-common (= 2.0.1-2+precise1) but it is
 not going to be installed
           Depends: libqgis2.0.1 but it is not going to be installed
           Depends: sip-api-8.1 but it is not installable
 qgis : Depends: libgdal1-1.7.0 but it is not installable
        Depends: libqgis2.0.1 but it is not going to be installed
        Depends: qgis-providers (= 2.0.1-2+precise1) but it is not going
        to be installed
        Recommends: qgis-plugin-globe but it is not going to be installed
 qgis-plugin-grass : Depends: libgdal1-1.7.0 but it is not  installable
     Depends: libqgis2.0.1 but it is not going to be installed
     Depends: grass641 but it is not installable
 E: Unable to correct problems, you have held broken packages.

Not sure what to do on this one


texts: Remove unwanted outline on raster export

texts:
I'm having an issue whereby the SVG I export via print composer doesn't reflect what I see in QGIS. My exported SVG image collects additional grey borders around the outermost boundaries of the raster.

UPDATE

Here I'm having the same problem with a different raster, on a different computer. This is QGIS Release 2.18. You can see my workflow below.

Exporting to SVG or PDF using print composer has a grey border each time


texts: Creating a Collector map that allows new records to be added to a related table while not allowing spatial features to be created or deleted?

texts:
I am trying to create inspection maps for Collector that are very similar to what can be found at this link:

https://blogs.esri.com/esri/arcgis/2015/02/18/related-tables-exploring-new-ways-to-use-collector-for-arcgis/

In order for these maps to function correctly for our needs, I need to have editing capabilities that would only allow for several attributes in a spatial layer to be edited and be able to add new inspection records to a related table (defined by a relationship class).  Except for this to work appropriately I cannot have the spatial layer so that it can have features added or deleted.

Is there a way to have this type of functionality while using Collector for ArcGIS?


texts: CORINE LandCover: what about the accuracy?

texts: How accurate is CLC? Do anyone know any scientific papers about that question? I tested the new CLC2006 for a region in turkey and its totally rubbish!


texts:
I am adding wms layers to my google maps application. As Google Maps projection is a web mercator, I thought it would be better to request wms layers with srs=3857 (or equivalent) instead of the 4326. I searched a bit for fellow implementation and I saw that most (if not all) use 4326 instead.

I wanted to experiment by myself, I request the 3857 to the map server and do some reprojections for the BBOX coordinates.
I have the two running so I can compare them. The result is:
- 4326, the tiles look OK acvtually, I can't notice much of a shift.
- 3857 the wms tiles are completely off, huge lat shift

My question is: I am right to think that I should be querying for a 3857 projected wms layer or should I stick to 4326? and is it normal that I have a big shift in my 3857 tiles?

edit

This is the code I use for creating the layer

    layer = new google.maps.ImageMapType({
        getTileUrl: function (coord, zoom) {
            var proj = _map.getProjection();
            var zfactor = Math.pow(2, zoom);
            // get tile coordinates
            var top = proj.fromPointToLatLng(new google.maps.Point(coord.x * 256 / zfactor, (coord.y + 1) * 256 / zfactor));
            var bot = proj.fromPointToLatLng(new google.maps.Point((coord.x + 1) * 256 / zfactor, coord.y * 256 / zfactor));
            // convert from wgs84 lat,lng to new srs coordinates
            var topConverted = srsConversion(top);
            var botConverted = srsConversion(bot);
            //create the Bounding box string
            // handles 1.3.0 wms by ordering lat, lng instead of lng, lat
            var bbox;
            if(version == "1.3.0") {
                bbox = botConverted.lng + "," + botConverted.lat + "," + topConverted.lng + "," + topConverted.lat;
            } else {
                bbox = topConverted.lng + "," + topConverted.lat + "," + botConverted.lng + "," + botConverted.lat;
            }
            //base WMS URL
            var url = lURL + "&BBOX=" + bbox; // set bounding box
            return url; // return URL for the tile
        },
        isPng: true,
        tileSize: new google.maps.Size(256, 256)
    });


texts: Line of Sight between sets of two points each

texts:
So I have a list of points with their longitude and latitude (Lon_A, Lat_A and Lon_B, Lat_B). I also have a DEM layer from the SRTM for the relevant location. I need to check if line of sight exists between these point pairs. What I'm essentially doing is checking if it's possible to create a microwave link between given point pairs.

So in a nutshell If I have a data set of point pairs I need to get a Boolean output whether Line of Sight exists or not. Wondering if there's a way to do this using QGIS or GRASS GIS. Looking for a process where this can be done en masse instead of checking each individual point pair.


texts: Scaling up an image and summing band values (rather than resampling)

texts:
I am trying to take an image at a native resolution scale of 250m and reproject it to 5000m and sum the values in the smaller pixels to the larger one.

I can do this in Earth Engine with the following, but it's just re-sampling the data (I'm guessing it's just calculating an average value):

var imageForAnalysis = ee.Image('JRC/GHSL/P2016/POP_GPW_GLOBE_V1/2015');
var proj = imageForAnalysis.projection();

var reprojectedImage = imageForAnalysis
    .reduceResolution({
      reducer: ee.Reducer.sum(),
      maxPixels: 1024
    })
    .reproject({
      crs: proj,
      scale:5000
    });

I've tried changing the reducer to ee.Reducer.sum(), but it does not yield a total value for the larger pixel size.

Any ideas?


texts: I have downloaded ALOS-1 PALSAR level 1.0 data from data archive  from Alaska Satellite Facility. After unzipping files ,I am getting image files,Leader Files,Trailer File and Volume File. All the files have extension of .0__A. I am trying to open the Image file(IMG) through gdal.Open and ReadAsArray() function but gdal is not able to open the data. Kindly suggest what can be done to visualise the data.I am using gdal 2.0.0.


texts:
I'm having trouble with the dissolve function in QGIS. Is it because of an overlap? It's not a GRASS layer, just shapefiles. I'm new with GIS software in general. Searched for similar problems but didn't find anything "monkey proof" so the speak.

Is it because of overlapping polygons? Geometry validity shows multiple errors but can't use v.clean because it's not a GRASS layer.


texts: How do you Select an existing vector Layer and Edit it, using python?

texts:
I'm developing a small plugin that needs to be able to select a vector layer programatically based on what layer the user has selected. However I'm not sure how to do so.

I believe I found the two elements I need to make it work: setCurrentLayer and selectedLayer from the QgisLegendInterface, but I'm not certain how they work or how they are implemented, I cannot either find any examples of them other than in the API.

Am I close? Or is there a better/proper form of doing this?


texts: coordinates of objects from geo-tagged images

texts:
I have geo-tagged images (i.e have the latitude and longitude from where the images were taken) from a truck every 20 feet.

I am trying to understand if it is possible to extract the latitude and longitude of the objects seen in the images.


texts: Sources for Free High-Resolution (4m or below) Orthophotography (Aerial Photography)

texts:
Does anyone know of reliable sites that provide high resolution Orthophotography for The USA?

I am not looking for Satellite Imagery. This is limited to Aerial Photography (Orthophotography) or Ortho Drone imagery in raster format - The data will be for Satellite Imagery Orthorectification.


texts: Serving 1 band GeoTIFF float32 file in Titiler/riotiler

texts:
I am very new to GIS and setting up COG servers.
The problem I am facing:
I am unable to serve my DEMs/DSMs, which are in GeoTIFF format, as .png map tiles (z/x/y.png) using a tiling server "TiTiler", which is based on riotiler.
However for normal RGB raster TIFFs, it works fine. I am not at all sure how to interpret DEM/DSM COGs to tiling server.
Titiler has an API endpoint which returns the GeoTIFF information and type. See this is the output when I feed in both RGB TIFF and my DEM:
DEM:
{
    "bounds": [88.42379642445006, 22.56651317616363, 88.44287051283406, 22.58178886547563],
    "band_metadata": [
        [1, {}]
    ],
    "band_descriptions": [
        [1, "band1"]
    ],
    "dtype": "float32",
    "colorinterp": ["gray"],
    "nodata_type": "Nodata"
}

RGB TIFF:
{
    "bounds": [88.42379642445006, 22.56651317616363, 88.44287045398606, 22.58178886547563],
    "band_metadata": [
        [1, {}],
        [2, {}],
        [3, {}],
        [4, {}]
    ],
    "band_descriptions": [
        [1, "band1"],
        [2, "band2"],
        [3, "band3"],
        [4, "band4"]
    ],
    "dtype": "uint8",
    "colorinterp": ["red", "green", "blue", "alpha"],
    "nodata_type": "Alpha"
}

Clearly the DEM file is a single band TIFF, and also float32 image file. Now the issue I am facing here is, Titiler (map tile server), is not able to generate PNG files based on float32 TIFF.
As of now I am really confused how to go about this. My only goal is to generate .png maptiles from my DEMs/DSMs, so that I can overlay them on Google Maps.


texts: Converting 32 bit DEM to 16 bit in ArcGIS Desktop?

texts: I have a 32 bit DEM with resolution 0.5 * 0.5 and range of elevation is 1299 to 1489. When I want to export or use copy raster tool or create fill in Arcgis, the results are elevation between -3.40282e+038 and 3.4282e+038 and the map does not show any elevation and range of color. I tried to convert 32 bit to 16 bit but it didn't work.


texts:
I have two versions of Python gdal running in different environments. One is slightly older (2.1.1) and gdal.Warp seems to work just fine there, the newer version doesn't. Has anyone else had this problem?

Edit:
I'm calling gdal.Warp like this:

gdal.Warp(output_raster, input_raster, dstSRS='EPSG:3083')

Error traceback, in the env with the newer version, looks like this:

runfile('C:/Users/mrt2594/clip&reproj_GRACE_TX.py', wdir='C:/Users/mrt2594')
Traceback (most recent call last):

  File "<ipython-input-7-b1b650156b3c>", line 1, in <module>
    runfile('C:/Users/mrt2594/clip&reproj_GRACE_TX.py', wdir='C:/Users/mrt2594')

  File "C:\Anaconda\envs\osgeo-lab\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 786, in runfile
    execfile(filename, namespace)

  File "C:\Anaconda\envs\osgeo-lab\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 110, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "C:/Users/mrt2594/clip&reproj_GRACE_TX.py", line 33, in <module>
    gdal.Warp(out_savepath, src_tif, options = clip_specs)

  File "C:\Anaconda\envs\osgeo-lab\lib\site-packages\osgeo\gdal.py", line 579, in Warp
    return wrapper_GDALWarpDestName(destNameOrDestDS, srcDSTab, opts, callback, callback_data)

  File "C:\Anaconda\envs\osgeo-lab\lib\site-packages\osgeo\gdal.py", line 3189, in wrapper_GDALWarpDestName
    return _gdal.wrapper_GDALWarpDestName(*args)

TypeError: in method 'wrapper_GDALWarpDestName', argument 4 of type 'GDALWarpAppOptions *'


texts:
I am fairly new to python scripting.
I am attempting to do some batch conversions from ASCII to Raster in order to save myself some time. However I am continually getting issues with my script.

My current script is:

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

# set current workspace
env.workpace = "G:/Arc_GIS/ANUCLIM_DATA/2020-2039/c-ECHAM5/R1/extract"
outworkpace = "G:/Arc_GIS/ANUCLIM_DATA/2020-2039/c-ECHAM5/R1/Python_out"
dataType = "INTEGER"

#Check out the arcGIS spatial Anaylst extension line
arcpy.CheckOutExtension("Spatial")

#Get a list of asciis
for file in arcpy.ListFiles("*.asc"):

    #Print raster list
    print file
    outputraster = file[1:36] + ".tif"
    #save TIF
    #outraster.save(output)
    arcpy.ASCIIToRaster_conversion(file,outputraster,dataType)

I am returning with the error

Traceback (most recent call last):
  File "G:/Arc_GIS/Coding/ASCIIToRaster_conversion", line 15, in <module>
    for file in arcpy.ListFiles("*.asc"):
TypeError: 'NoneType' object is not iterable

Where am I going wrong?


texts: Deleting floating lines from network using GRASS 7 in QGIS

texts:
I'm trying to find a way to remove floating line segments from a road network shapefile.

Is there a relatively fast way to identify and delete these features?

I've been using GRASS 7 in QGIS to remove dangles < 100m and generalize its shape, but I haven't found an appropriate tool for removing these segments that are not part of the main network.

I'm open to finding solutions in QGIS, but I'm curious to see what GRASS tools can do for situations like this.


texts: Why pgr_dijkstra() only finds the routes that are pretty straight?

texts:
I am trying to use pgRouting in order to find the shortest route between two points. I realized that pgr_dijkstra() function is not able to find the route when  the two points are far from each other or when they are not very straight. However, when I select the two points in a neighbourhood, it can give me the shortest path.

This is the roads table (edges):

and this is the table of nodes that pgrouting created with the ending of "_vertices_pgr":

Does any one know what is the problem?
I used this query:

select pgr_dijkstra('select id , source, target, st_length(geom) as cost FROM q_roads_noded', 10237, 10551)

q_roads is the table of roads (edges)

The network is like this: I tried to find the shortest path between points with id: 10237 and 10551. But the function returned nothing!!!


texts: How to programmatically remove a donut hole from Polygon in arcobjects


texts: how to batch reclassify rasters with percents range

texts:
Using Model Builder I need to reclassify slope rasters with percents range value. I have read and ask here batch reclassify rasters with different elevation ranges, but I still need more explanation.

Each one of my rasters has different ranges, this is an example from manual classify table:

**Image 1**
Old Values,New Values
0-278,1 (8%)
278-523,2 (15%)
532-871,3 (25%)
871-1394,4 (40%)
1394-3486,5 (100%)

**Image 2**
Old Values,New Values
0-471,1 (8%)
471-883,2 (15%)
883-1472,3 (25%)
1472-2355,4 (40%)
2355-5889,5 (100%)

I need to keep that 8%,15%,25%,40% & 100% ranges. I'am thinking to read each image ranges and write it to a table on the fly, but I don't know for sure how to do that.

Edit:

I think I find a way to do it, just not sure how it is done.

1. I'll get the maximum value using Get Raster Properties tool
2. then do a mathematic calculations (maxValue "multiplied by" (8 "divided by" 100))
3. write the output to dbf file with a format Low,High,Val (fill in 0 as the first Low field)
4. iterate the mathematic calculation for the next percents range (15%,25%,40% and 100%)
5. Fill in the last result of the mathematic calculation in the Low field of an on going calculation (example below)
6. load the dbf file to the Reclassify by Table tools, DONE.

I still need help with mathematic calculation part, I don't know Python.

Example:
MaxValue = 100

Low, High, Val
0, 8, 1
8, 15, 2
15, 25, 3
25, 40, 4
40, 100, 5


texts: Inserting image in an attribute table in QGIS Composer?

texts:
I'm trying to create second page of map atlas. In this page I need to show all records in the range of current atlas page.

I've used Insert Attribute Table tool in composer and it worked. But I need to add a graphic to every record (on the end on attribute table).

Is there a possibility to show image from attribute table in composer or I need to use some separate HTML code?


texts: SQL statement to find differences between two fields

texts: I am trying to find the difference of 11 feet or more between two fields (Both are Double types). For example, I have a field called shape_length and another field called station_length and I want to see if there are any differences of 11 feet or more between the two, and select those differences. I think it's a simple SQL statement in the "Select By Attributes dialog in ArcGIS but I can't figure it out.


texts: Open Image from Sentinel 1 in SNAP problems

texts:
I have download SNAP today but I'm facing few problems with upload the image.
I have download image from Copernicus Open Acess Hub but whenever I download it I can't open it properly.
those are the filles that were downloaded-

in measruemnet there are two images -vv and vh, as I download, but  as two different images.
When I open one of it in SNAP, I get one image with one band that looks like this-

whenever I try to run any process on the image I get this error message-

how can I solve this and upload image to the SNAP?


texts: FME conversion defaults to inches when all data is in meters

texts:
I regurlaly do conversions between shape, dwg, tif, and a few other formats in FME. Often, the resulting unit in dwg files is inches. I'm a bit surprised at this, since all ingoing data is in meters and also in a coordinate system defined in meters.

Is there a way to set the units in FME to make sure it's using meters in the dwg file?

This is a recent workspace I have used to create a dwg file. All data is set to Sweref99_1800 (EPSG3011), which is based on meters. The dwg file was in inches when done.

I'm guessing this is a setting for default values in Autocad, but I'm hoping I can access the setting from FME. I'm also open to alternative solutions, in case I have got something wrong! :)


texts: Can we use "evaluate" function for other models not included in dismo package?

texts:
I'm using several models to predict species occurence probabilities. Some of them are not included in dismo library such as "MARS" and random forests. I want to know if I still can use "evaluate" function of dismo to evaluate those models.  I mean when I run this line :
eval_mars_test <- evaluate(p=envtest_p , a=envtest_a,  model= model_mars)

It shows no error but how can I confirm that TPR and FPR here are correctly calculated?


texts: How can I create individual directional point symbols in MapInfo?

texts: I have a point tab file where one of the attributes is a 0-360 direction. I want to create an arrow symbol for  each point with the arrow pointing in the direction listed for that point in the table. Is this possible? I'm using MapInfo but would welcome a solution in ArcMap or QGIS also.


texts: python psycopg2 insert table into postgis: geometry requires more points

texts:
I have a feature class of 20k records. trying to import it into postgis using arcpy and psycogp2. I have a pretty basic script

cur.execute('''drop table if exists developed_lands.daysemtric_test''')
cur.execute('''CREATE TABLE developed_lands.daysemtric_test
            (pid serial primary key,
            pop_dns float,
             shape geometry);''')
with arcpy.da.SearchCursor(shp, ["POP_DNS","SHAPE@WKT"]) as cur2:
    for row in cur2:
        sql = "insert into developed_lands.daysemtric_test(pop_dns,shape)values(%s,ST_GeomFromText(%s)"
        cur.execute(sql,(row[0],row[1]))
        x+=1
        print x
conn.commit()

error:

   Traceback (most recent call last):
  File "C:\Users\rizagha\Desktop\module1.py", line 16, in <module>
    cur.execute(sql,(row[0],row[1]))
InternalError: geometry requires more points
HINT:  "...6391.65285167098 737288.47085659206)," <-- parse error at position 13049795 within geometry

                                                                   ^

This script works because If I skip the first feature it inserts the records up to row 2,000 then throws the same error. I put the function st_makevalid() around ST_GeomFromText(%s) and it did nothing to help. running repair geometry on the feature class now, but that usually takes hours...

the first row where the error is happening, the geometry is extremely complex. is it possible the WkT is too long?

UPDATE

repair geometry finished and I reran the script and I got the same error