texts: Styling polygons in a map tile PBF file in Leaflet

texts:
How can I dynamically add data to a pbf file generated by Maptiler? I created a pbf file with the MapTiler desktop application. I started with a geojson file containing polygon geometries. My file is stored on Maptiler cloud. Maptiler cloud provides this URL to include in my client code.

"https://api.maptiler.com/tiles/28815aac-ebdc-4442-a4b3-62172c26ca52/{z}/{x}/{y}.pbf?key=<myMaptilerKey>";

I am importing the file into my leaflet map like this...

const FOO_PBF =
  "https://api.maptiler.com/tiles/28815aac-ebdc-4442-a4b3-62172c26ca52/{z}/{x}/{y}.pbf?key=<myMaptilerKey>";
  var vectorTileOptions = {
    interactive: true
  };
  var fooLayer = L.vectorGrid.protobuf(
    FOO_PBF,
    vectorTileOptions
  );

  var overlays = {
    FooLayer: fooLayer
  };

The polygons appear correctly. However, I would also like to add data to this layer. Specifically, I have a statistic pertaining to each polygon. I'd like to color the polygons with colors that represent the value of that statistic.

In the past, I imported the geojson file directly into my leaflet map and achieved this "chloropleth" type styling like this.

  var FooLayer= new L.geoJson.ajax(
    "https://raw.githubusercontent.com/...path_to_geojson_file_geo.min.json",
    {
      // style: defaultPolygonStyle,
      onEachFeature: function(feature, layer) {
        layer.bindPopup(
          feature.properties.statisticValue +
            ", " +
            $this.zipData[feature.properties.statisticValue]
        );
        layer.setStyle({
          color: getColorFromRedBlueRange(
            $this.zipData[feature.properties.statisticValue]
          ),
          weight: 0,
          fillOpacity: 0.3
        });
      }
    }
  );

where

  function getColorFromRedBlueRange(d) {
    return d > 0.8
      ? "#FF0000"
      : d > 0.7
      ? "#FF006F"
      : d > 0.55
      ? "#FF00EF"
      : d > 0.45
      ? "#DE00FF"
      : d > 0.3
      ? "#BC00FF"
      : d > 0.2
      ? "#6600FF"
      : "#00009FF";
  }

How can I similary style the polygons in my pbf file? The reason we are using the pbf file rather than the raw geojson is to improve performance by reducing network payload sizes.

With Stefan's tips, I found that this applies the specified yellow fill and red outlines to the polygons in the pbf tile layer...

var vectorTileOptions = {
        interactive: true,
        // pane: "OverlayPane",
        vectorTileLayerStyles: {
          layer_name: {
            fillColor: "yellow",
            fill: true,
            color: "red"
         }
        }
      };

However, if I try to style each polygon according to a statistic in the pbf file, this does not work...

  function getColorFromRedBlueRange2(d) {
    return d == 8
      ? { fillColor: "#FF0000", fill: true, opacity: 1 }
      : d == 7
      ? { fillColor: "#FF006F", fill: true, opacity: 1 }
      : d == 6
      ? { fillColor: "#FF00EF", fill: true, opacity: 1 }
      : d == 5
      ? { fillColor: "#DE00FF", fill: true, opacity: 1 }
      : d == 4
      ? { fillColor: "#BC00FF", fill: true, opacity: 1 }
      : d == 3
      ? { fillColor: "#6600FF", fill: true, opacity: 1 }
      : { fillColor: "#0009FF", fill: true, opacity: 1 };
  }
  var vectorTileOptions = {
    interactive: true,
    // pane: "OverlayPane",
    vectorTileLayerStyles: {
      layer_name: properties => {
        getColorFromRedBlueRange2(properties.SLDUST % 10);
      }
    }
  };

Why not? The SLDUST variable does exist in the properties object and a console.log statement does confirm that this code pulls the SLDUST values from the pbf file.

In both cases, I am adding the layer to the map with this code chunk...

  var pbfLayer = L.vectorGrid.protobuf(
    FOO_PBF,
    vectorTileOptions
  );


texts: How to update OSM data, without recreating the PostGIS database?

texts:
We are currently working with GeoServer and PostGis which stores the OSM data for Switzerland.

I would like to know, how to update the OSM data for my GeoServer. At my current level of knowledge, I would need to download the diff file than update the pbf file with osmosis after that I would have to create a new database for the new data and import it again.

With this method, I would need to reaply all configuration on my GeoServer again from the beginning.

Is there an easier way to update my OSM data?


texts: Altering label styles (font, color, size, etc.) using ArcPy?

texts:
I am trying to write a python macro for ArcGIS 10.2 that will automate the process of labeling points in a layer based on different types of apartments.  For example, one style of apartments would have a blue square with white bold text as one label; another would have a pink square with black bold text; another would have a red square with a white bold text, etc.

What expression or function would you use in arcpy to alter the style of these labels?

Currently, the script labels the apartments by a map code and changes its name in the TOC to apartments, as it should.

Now I want to change the font, graphic, and size of the labels to look make a unique label for each, and then make a button to press and simply automate the process any time we have to do one of these layers.

Is there any way the fonts, labels, and graphics can be changed with the use of Python scripting (ArcPy, IDLE, etc.)?

I decided to edit this post to give a final description of what I'm going for here so there's no confusion.  Essentially, the points on the map will go from looking like the above picture to this:

I have tried using the lblClass.expression function, however the code I am trying to use doesn't work.  I figure that is probably the step in the right direction. Below is some code that I've tried to use just to get the labelClass.expression to change the fonts of the labels.

for lblClass in aptLayer.labelClasses:
lblClass.showClassLabels = True
def getLabel (Code, Type):
if Type == "MRR":
    return "<FNT name = 'Arial' size = '12'>" & Code & "</FNT>"
elif Type == "TAX":
    return "<FNT name = 'Arial' size = '12'>" & Code & "</FNT>"
elif Type == "MRT":
    return "<FNT name = 'Arial' size = '12'>" & Code & "</FNT>"
elif Type == "GSS":
    return "<FNT name = 'Arial' size = '12'>" & Code & "</FNT>"
elif Type == "TGS":
    return "<FNT name = 'Arial' size = '12'>" & Code & "</FNT>"
else:
    return "<FNT name = 'Arial' size = '12'>" & Code & "</FNT>"

lblClass.expression = getLabel( [Map_Code],[Type])


texts:
I am trying to use Python in ArcGIS Pro to cleanse some data via imported Excel sheets. The expression works, but the Code Block if-else statement to return nothing if !Street_Id! is  returns an error.
Address
str(!Street_ID!) + ' ' + !Street!.split("/")[0]

Code Block
def f(x): ' '
if !Street_ID! = 'Null'
return x
else return str(!Street_ID!) + ' ' + !Street!.split("/")[0]

Error
File"<string>",line 2
if!Street_ID! = 'Null'
   ^
SynatxError: invalid syntax


texts:
I am trying to show the labels on the map on one particular field (Candidate Code) and its a String value, and also trying to set the scale range and placement property on labels (show labels out beyond 1:1500000 and top center only) But I am not sure how to do it.

Below is my complete code:

import arcpy
from arcpy import env
import arcpy.mapping

# Set environment settings
env.workspace = "C:\Users\Documents\Regions"
arcpy.env.overwriteOutput = True

#define projection
Prjfile = r"C:\Users\Documents\Regions\WGS1984.prj"
spRef = arcpy.SpatialReference(Prjfile)

#Local Variables
FC6 = r"C:\Users\Documents\Regions\West_DO.gdb\Candidate_Fbr_500M"
FC7 = r"C:\Users\Documents\Regions\West_DO.gdb\Search_Ring"

mxd = arcpy.mapping.MapDocument(r"C:\\Users\\Documents\\Demo_West.mxd")
df = arcpy.mapping.ListDataFrames(mxd,"Layers")[0]

arcpy.MakeFeatureLayer_management(FC7, r"C:\Users\Documents\Regions\Search_Ring.lyr")
layer7 = arcpy.mapping.Layer(r"C:\Users\Documents\Regions\Search_Ring.lyr")
Sym_layer7 = r"C:\Users\Documents\Symbology layer new\Search Ring.lyr"
arcpy.ApplySymbologyFromLayer_management (layer7, Sym_layer7)
arcpy.mapping.AddLayer(df, layer7, "AUTO_ARRANGE")
lyr7 = arcpy.mapping.ListLayers(mxd, "", df)[0]
lyr7.visible = False
lyr7.name = "Search_Ring"
lyr7.maxScale = 0;
lyr7.minScale = 2000000;
lyr7.transparency = 25;

def FindLabel ( [Candidate_Code] ):
    return  "<CLR blue='255'><FNT size = '9'><BOL>" + [Candidate_Code] + "</BOL></FNT></CLR>"

lyr7.showLabels = True

arcpy.MakeFeatureLayer_management(FC6, r"C:\Users\Documents\Regions\Candidate_Fbr_500M.lyr")
layer6 = arcpy.mapping.Layer(r"C:\Users\Documents\Regions\Candidate_Fbr_500M.lyr")
Sym_layer6 = r"C:\Users\Documents\Symbology layer new\Candidate_Fbr_500M.lyr"
arcpy.ApplySymbologyFromLayer_management (layer6, Sym_layer6)
arcpy.mapping.AddLayer(df, layer6, "AUTO_ARRANGE")
lyr6 = arcpy.mapping.ListLayers(mxd, "", df)[0]
lyr6.visible = False
lyr6.name = "Candidate_Fbr_500M"

mxd.saveACopy(r"C:\\Users\\Documents\\West.mxd")
print "save mxd"

arcpy.RefreshActiveView()
arcpy.RefreshTOC()
del mxd,df


texts:
I have 3 point feature classes (test30, test90, test180) in a scratch.gdb. I would like to create a kernel density raster for each dataset. I would like the output of the kernel density raster to be saved in the scratch.gdb.

The error I am receiving is ERROR 010093: Output raster format UNKNOWN is unsupported. Can I save a kernel density raster to a File Geodatabase?

I am new to python and functions, so I may have some of this logic written incorrectly.

ArcDesktop 10.3.1

import os, sys, re, shutil
import arcpy
from arcpy.sa import *
import arceditor

cwd = os.getcwd()

scratchGDB = r'C:\Users\test\scratch.gdb'
arcpy.env.workspace = scratchGDB
arcpy.env.overwriteOutput = True

# kenrnel density
def kernel_density(dataset):

    try:
        if arcpy.CheckExtension("Spatial") == "Available":
            arcpy.CheckOutExtension("Spatial")
        else:
            raise LicenseError

        kd = KernelDensity(dataset, population_field='NONE', cell_size=1000, search_radius=50000, area_unit_scale_factor='SQUARE_MILES')

        kd.save(scratchGDB)

    finally:
            arcpy.CheckInExtension("Spatial")

kernel_density('test180')
kernel_density('test90')
kernel_density('test30')


texts: How to resample a DEM from GCS to UTM with respect to reprojection bounds parameter?

texts:
When reprojecting a DEM (e.g., SRTM-derived 1-arc DEM) using PCI Geomatica for example, to constrain pixel size and maintain a square pixel you can either change the image extent  or the number of pixels, or both.
What are the implications of changing either parameters with respect to preserving the DEM accuracy?


texts: Can I specify a point and get all features within a certain radius of this point?

texts:
I have recently shifted to QGIS from ArcGIS and am looking for some help/pointers with the following:

Does QGIS have a method by which I can specify a lat-lon coordinate (or a utm x,y coordinates, coordinate system not imp), and get all features that lie within a certain radius of this point?

I downloaded the 3rd party plugins for QGIS and after looking through all of them, didn't find anything with this functionality (though I could have very well not understood something completely). I am considering writing a plugin in case something doesn't already exist

Edit: I looked at the option under: View -> Select -> Select Features By Radius. This seems to be similar in concept to what I am aiming to do (except that the center needs to be specified).


texts: Selecting specific row in attribute table using ArcPy (without iterating through whole table)?

texts:
I am using ArcGIS 10.2 and have assumed until now that it is possible to select a specific row in an attribute table. However, I can't seem to find a way to do this. The only available methods seem to be to return a search or update cursor, and to iterate through each row in the table.

Is there really no way to select only the desired row?

What I am trying to do is to select a target row - e.g. 50423 - and then to update a field for that row. I'd hate to iterate through the entire table of millions of records.


texts: How to extract the node and cell data values in QGIS?

texts:
I have just started using GIS applications so i am a novice in this field. I am currently using QGIS. The project demands the unstructured(triangular) meshing which i was able to achieve through QGIS-GMSH plugin. Now I have it meshed and next I want is to extract the node(x, y coordinates of all nodes) and cell data(nodes comprising the triangular cell and its coordinates) values and that's where i am stuck. I googled it but no luck.

I hope I have used the right terminology and pardon me for any mistakes.


texts: Convert X,Y State Plane coordinates to decimal degrees

texts: I have 1000's of data point I need to convert from North Central Texas 4202 into Lat/Long decimal degrees. Can someone give me a conversion equation so I can quickly convert using excel? I have only been able to find online conversion calculator which convert one entry at a time. This is not time efficient for me. Please help!


texts:
I have multiple WMS in OpenLayers with a layer switcher and can be turned on and off. Now I have an image of each WMS, what I want is to show image of each WMS on click OF that corresponding WMS. Every WMS has only one polygon. How that can be done?

I apologize if it is a very basic question.

Below is my JavaScript code:

function initMap()
  {
  var osmLayer = new OpenLayers.Layer.OSM();
  var gmap     = new OpenLayers.Layer.Google("Google Streets", {visibility: false});
  var info;
var Guyana= new OpenLayers.Layer.WMS("Guyana Country",
 "http://localhost:8080/geoserver/wms",
  {
    format: "image/png8",
    transparent: true,
    layers: "opengeo:guy_adm0",

  }, {
    opacity: 0.6,
    singleTile: true,
    visibility: true,
  });
var Chile= new OpenLayers.Layer.WMS("Chile Country",
 "http://localhost:8080/geoserver/wms",
  {
    format: "image/png8",
    transparent: true,
    layers: "opengeo:chl_adm0",

  }, {
    opacity: 0.6,
    singleTile: true,
    visibility: true,
  });
var Venezuela= new OpenLayers.Layer.WMS("Venezuela Country",
 "http://localhost:8080/geoserver/wms",
  {
    format: "image/png8",
    transparent: true,
    layers: "opengeo:ven_adm0",

  }, {
    opacity: 0.6,
    singleTile: true,
    visibility: true,
  });
var Columbia= new OpenLayers.Layer.WMS("Columbia Country",
 "http://localhost:8080/geoserver/wms",
  {
    format: "image/png8",
    transparent: true,
    layers: "opengeo:col_adm0",

  }, {
    opacity: 0.6,
    singleTile: true,
    visibility: true,
  });
var Argentina= new OpenLayers.Layer.WMS("Argentina Country",
 "http://localhost:8080/geoserver/wms",
  {
    format: "image/png8",
    transparent: true,
    layers: "opengeo:arg_adm0",

  }, {
    opacity: 0.6,
    singleTile: true,
    visibility: true,
  });
var Bolivia= new OpenLayers.Layer.WMS("Bolivia Country",
 "http://localhost:8080/geoserver/wms",
  {
    format: "image/png8",
    transparent: true,
    layers: "opengeo:bol_adm0",

  }, {
    opacity: 0.6,
    singleTile: true,
    visibility: true,
  });
var Paragua= new OpenLayers.Layer.WMS("Paragua Country",
 "http://localhost:8080/geoserver/wms",
  {
    format: "image/png8",
    transparent: true,
    layers: "opengeo:pry_adm0",

  }, {
    opacity: 0.6,
    singleTile: true,
    visibility: true,
  });
var myWord = "SNOW COVER BASINS NORTHERN AREAS OF PAKISTAN";

// Map with projection into (required when mixing base map with WMS)
  olMap = new OpenLayers.Map("divid", {
  projection: "EPSG:900913",
   displayProjection: "EPSG:32463",
  units: "m",
  layers: [gmap,osmLayer,Guyana,Chile,Venezuela,Columbia,Argentina,Bolivia,Paragua]
});

olMap.addControl(new OpenLayers.Control.LayerSwitcher());
olMap.setCenter(new OpenLayers.LonLat(6011274.179826, -5256566558465465.6),3);

    }


texts: Analytical hierarchy process tool - SAGA. How to configure Pairwise comparisons table?

texts: What is the configuration and extension (.txt, .csv ...) of the Pairwise comparisons table file requested by the SAGA AHP tool- Analytical Hierarchy Process - in QGIS? Are there examples of applications available?


texts: How can I display annotations for certain features?

texts:
I am trying to display annotations for single features. I know it is possible to display annotations for layers, and for that entire map, but I would like to display them based on a SINGLE selected feature.
example:-

I have a layer called A that is visible on the map, so all the features in it are visible. Now I select a feature, and want to display the annotation for that particular feature not the whole layer; how is this possible?

In my code I am doing this but for the entire layer, but I want it for a single feature, any ideas?

Here is my code which is a method that displays annotations for STRUCTURE features:

public static void ToggleXyLabelVisibility(IApplication application,ISchematicLayer layer, ISchematicFeature slctdFtr =null) {
        var doc = application.Document as IMxDocument;
        if (doc != null) {
      // display annotation for the entire layer if no features are selected
            if (slctdFtr == null) {
                var compLayer = layer as ICompositeLayer;
                if (compLayer != null) {
                    for (int i = 0; i < compLayer.Count; i++) {
                        var geoflay = compLayer.Layer[i] as IGeoFeatureLayer;
                        if (geoflay != null) {
                            for (int j = 0; j < geoflay.AnnotationProperties.Count; j++) {
                                var annotationProperties = geoflay.AnnotationProperties;
                                IElementCollection col1;
                                IElementCollection col2;
                                IAnnotateLayerProperties properties;
                                annotationProperties.QueryItem(j, out properties, out col1, out col2);
                                if (string.Equals(properties.Class, Locationlabel,
                                    StringComparison.CurrentCultureIgnoreCase)) {
                                    properties.DisplayAnnotation = !properties.DisplayAnnotation;

                                }
                            }
                        }
                    }
                }
             }
            // if features are selected, then display it for those selected features only
            else {
            // CODE GOES HERE WHICH I CAN'T FIGURE OUT HOW TO DO
            }
            doc.ActiveView.Refresh();
        }
    }


texts: Cannot check in related tables in ArcPad?

texts:
I'm using ArcPad Data Manager (ArcMap 10) to check in data in axf file. I can check in feature layers, but I can't check in related tables and the feature class parent of the related tables.

I have tried changing the [AXF_STATUS] to 1 in the parent feature layer and in its child related table, but it doesn't work too.


texts: Qgis: Export only some columns of the selected objects to CSV

texts:
I have a layer with a table containing several columns.
If I right click on the layer --> Save as.. --> CSV and chose to only include my selection I get a nice CSV file with all columns for the selected objects.

However, I only want the data from one specific column to be exported to the CSV file.

How would I go around this?


texts: How to export Std Dev from Raster Dataset of DEM to classify current extent and execute Contour Tool with Python?

texts:
I want to create Countour Lines for a DEM which are based on a "STD dev" classificaion.
My problem is that the "Std dev" is exported in an unknown type instead of a String. So the classifacation for the Contour Tool can not be executed.

Do you have any ideas how i could convert it in a String?

STD

elevstd=arcpy.GetRasterProperties_management("rectExtract","STD")

Define Variables

Intervall=0
Automatisch=0

If Statements

if Intervall==Automatisch and elevstd>"0" and elevstd<="5":   # Class1
    Contour("rectExtract", "C:/fakepath/Class1.shp", 1,0)
if Intervall==Automatisch and elevstd>"5" and elevstd<="10":  # Class2
    Contour("rectExtract", "C:/fakepath/Class2.shp", 5,0)

Thanks for your help,


texts:
I had developed a tool by using C# (ArcObjects) and I want to (import that customized DLL and) call a particular function from that class by using ArcPy.

How can I achieve this?


texts: Write information on file using exiftool_call in R

texts:
I have a series of images (.tiff) that are only georeference at the true centroid of the picture. I have calculated the missing coordinates (top right, top left, bottom right and bottom left) in the X and Y plain, and now I will like to write those coordinates in each image.

In order to do that, I was told that the exifr package in R has a function call exiftool_call that could be helpful to complete that task. I have tried different things, but they don't seem to be working, so I was wondering if anyone knows how to accomplish this task using exiftool_call

Here is what I have tried so far:

coordinates.data<-read.csv("...",header=TRUE)

list.pics <- list.files(path ="...",pattern = "*.tif",full.names = TRUE)

exiftool_call(args = coordinates.data, fnames = list.pics)

exiftool_call(args = c(coordinates.data$TopLeftx,coordinates.data$TopLefty,coordinates.data$TopRightx,
                      coordinates.data$TopRighty,coordinates.data$BottomLeftx,coordinates.data$BottomLefty,
                      coordinates.data$BottomRightx,coordinates.data$BottomRighty),fnames = list.pics)

exiftool_call(args = c("coordinates.data$TopLeftx","coordinates.data$TopLefty","coordinates.data$TopRightx",
                       "coordinates.data$TopRighty","coordinates.data$BottomLeftx","coordinates.data$BottomLefty",
                       "coordinates.data$BottomRightx","coordinates.data$BottomRighty"),fnames = list.pics)

Output str(coordinates.data)

'data.frame':   5 obs. of  16 variables:
 $ X                          : int  1 2 3 4 5
 $ meta.df.FileName           : Factor w/ 5 levels "IMG_0000_1.tif",..: 1 2 3 4 5
 $ meta.df.FileInodeChangeDate: Factor w/ 3 levels "2018:07:11 14:16:48-05:00",..: 1 3 1 2 2
 $ meta.df.BandName           : Factor w/ 5 levels "Blue","Green",..: 1 2 4 3 5
 $ coorMetric.lon             : num  329878 329878 329878 329878 329878
 $ coorMetric.lat             : num  4057023 4057023 4057023 4057023 4057023
 $ ImageWidth                 : num  88.8 88.8 88.8 88.8 88.8
 $ ImageHeight                : num  66.7 66.7 66.7 66.7 66.7
 $ TopLeftx                   : num  329833 329833 329833 329833 329833
 $ TopLefty                   : num  4057056 4057056 4057056 4057056 4057056
 $ TopRightx                  : num  329922 329922 329922 329922 329922
 $ TopRighty                  : num  4057056 4057056 4057056 4057056 4057056
 $ BottomLeftx                : num  329833 329833 329833 329833 329833
 $ BottomLefty                : num  4056989 4056989 4056989 4056989 4056989
 $ BottomRightx               : num  329922 329922 329922 329922 329922
 $ BottomRighty

   : num  4056989 4056989 4056989 4056989 4056989

One File Example

one.Image<-read_exif(path = "...")
one.Image

coordinate.one.Image<-read.csv("...", header = TRUE)
head(coordinate.one.Image)

exiftool_call(args = coordinate.one.Image, fnames = one.Image)

Output str(coordinate.one.Image)

'data.frame':   1 obs. of  15 variables:
 $ meta.df.FileName           : Factor w/ 1 level "IMG_0001_1.tif": 1
 $ meta.df.FileInodeChangeDate: Factor w/ 1 level "2018:07:11 14:16:49-05:00": 1
 $ meta.df.BandName           : Factor w/ 1 level "Blue": 1
 $ coorMetric.lon             : num 329860
 $ coorMetric.lat             : num 4056981
 $ ImageWidth                 : num 88.8
 $ ImageHeight                : num 66.7
 $ TopLeftx                   : num 329815
 $ TopLefty                   : num 4057015
 $ TopRightx                  : num 329904
 $ TopRighty                  : num 4057015
 $ BottomLeftx                : num 329815
 $ BottomLefty                : num 4056948
 $ BottomRightx               : num 329904
 $ BottomRighty               : num 4056948


texts: Spatial join with multiple intersecting features

texts:
I have a polygon dataset representing pavement areas broken up by year of construction and a roads centerline dataset that contains values for remaining service life (RSL).  I want to spatially join the two datasets so that the RSL values are added to the correct pavement polygons.  However the centerlines are not split along the same lines that the pavement polygons were so in some cases multiple lines intersect a single polygon.  To avoid this I though I would set the HAVE_THEIR_CENTER_IN match_option, in ArcMap 10.2.2, but upon further inspection there are shorter polygons that correspond to a longer line segment that will not have its center within the polygon.  This picture tells it better:

The arrows represent the breaks of the lines. Notice how some polygons have multiple lines intersecting them, and others will not contain the center point of the line segment.

My tentative plan is to split the line segments by the polygons and then remove the small line segments that do not correspond to the pavement polygon they fall into.  This will be tedious and I was hoping for a better way.  Open to QGIS, ArcMap and Python solutions.

Is there a way to do a spatial join where only the values from the longer line segments that intersect a polygon are kept while the shorter line segment(s) are ignored?


texts: How to determine if a polygon shapefile is multipart or singlepart?

texts: The ArcMap Properties dialog does not indicate if a shapefile contains multipart or singlepart polygons.  Is there some other way to identify if a polygonal shapefile (or any feature class, for that matter) contains singlepart or multipart?


texts: ImportError DLL load failed when import qgis.core but no errors when import qgis

texts:
EDIT
So.. after some problems everything seems to semi-work I guess ? I will update my variables so maybe it will help someone.

I am working on a standalone application using Python and QGIS.
I've read every single tutorial/forum I could find and tried multiple times the method described in the QGIS Python Programming CookBook but I can't make it work.
I am using :

Python 3.7.0
OSGeo4W64 3.12
PyQT5 5.15.1

I've already configure the PATH and PYTHONPATH variables using the commands below in the QGIS Python console (the path that I obtained where added to the variables)
import sys
print(sys.path)
qgis

Paths that didn't existed where ignored.
For the details here are my variables :
PYTHONPATH
 - C:\OSGeo4W64;
 - C:\OSGeo4W64\bin;
 - C:\OSGeo4W64\apps\qgis\python;
 - C:\OSGeo4W64\apps\qgis\python\qgis;
 - C:\OSGeo4W64\apps\Python37;
 - C:\OSGeo4W64\apps\Python37\lib;
 - C:\OSGeo4W64\apps\Python37\lib\site-packages;
 - C:\OSGeo4W64\apps\Python37\lib\site-packages\win32;
 - C:\OSGeo4W64\apps\Python37\lib\site-packages\win32\lib;
 - C:\OSGeo4W64\apps\Python37\lib\site-packages\pythonwin;
 - C:\OSGeo4W64\apps\Python37\DLLs;
 - C:\OSGeo4W64\apps\Python37\Scripts;
 - C:\Users\ [username]\AppData\Roaming\QGIS\QGIS3\profiles\default\python;

PATH
 - C:\OSGeo4W64\bin;
 - C:\OSGeo4W64\apps;
 - C:\OSGeo4W64\apps\qgis\bin;
 - C:\OSGeo4W64\apps\qgis\python;
 - C:\OSGeo4W64\apps\qgis\python\qgis;
 - C:\OSGeo4W64\apps\Python37;
 - C:\OSGeo4W64\apps\Python37\lib;
 - C:\OSGeo4W64\apps\Python37\lib\site-packages;
 - C:\OSGeo4W64\apps\Python37\lib\site-packages\win32;
 - C:\OSGeo4W64\apps\Python37\lib\site-packages\win32\lib;
 - C:\OSGeo4W64\apps\Python37\lib\site-packages\pythonwin;
 - C:\OSGeo4W64\apps\Python37\DLLs;
 - C:\OSGeo4W64\apps\Python37\Scripts;
 - C:\Users\ [username]\AppData\Roaming\QGIS\QGIS3\profiles\default\python;
 - C:\OSGeo4W64\apps\Qt5\bin;
 - C:\OSGeo4W64\apps\Qt5\lib;
 - C:\OSGeo4W64\apps\Qt5\plugins;

And here's the kicker :
>>> import qgis
>>>
>>> import qgis.core
Traceback (most recent call last):
   File "<stdin>", line 1 in <module>
   File "C:\OSGeo4W64\apps\qgis\python\qgis\core\__init__.py", line 25, in <module>
      from qgis._core import *
ImportError: DLL load failed: the specified procedure could not be found.
>>>

(I've translated the last sentence since it was in french)
So as you can see, when I simply import qgis, it works but if I import qgis.core / qgis.gui and so on, I have the same error.
I've read that changing the version of QGIS might help but I want to try something else before switching to another version and probably messing all up again.


texts: Turf Envelope returns 5 array coordinates for rectangle polygon

texts:
I'm trying to understand why turf.envelope() returns 5 coordinate arrays.  I'd expect to see only 4.  Is this normal/expected behavior?

I ask because I'm currently trying to debug an issue I'm having while setting a bounding box in MapBox after being supplied x amount of coordinate arrays from the draw control from the draw.update event and I'm wondering if this may be the culprit.

I'm leaning toward having 5 coordinate arrays being returned not being the issue, because everything functions as expected when processing the draw.create event.


texts: Displaying overlapping lines as separate offset lines in QGIS3

texts:
I'm drawing lines in a map that represent cables in underground pipes. Sometimes there are multiple cables going through the same pipe, so I'd like to draw each cable as a different line, but want it to be represented/rendered as separate lines having some offset and ideally even being differently colored. Have tried multiple approaches with symbology and geometry generator, but haven't really had a success. The main goal is to do it via rendering/display and not have to manually set different offset value to each line or draw parallel lines myself. Is this even possible in QGIS?
Let's say three separate cables are coming towards the manhole from the left, one cable goes up, two goes down:

I want them to be displayed something like this to visually represent their actual paths:


texts:
I am trying to oversample a part of a big GeoTIFF DEM to display it as a hillshade. I'm doing the process in 3 steps:

extract a part of my GeoTIFF
oversample that part
generate the hillshade with QGIS

Here's the command I use to extract a part:
gdal_translate \
    -co BIGTIFF=YES -co TILED=YES -co COMPRESS=LZW -co PREDICTOR=2 \
    -projwin 3.4999998866667914 43.12434922166657 5.968749886666799 40.78056666958323 \
    dem.tif \
    a.tif

When openning a.tif in QGIS and generating a hillshade (using bicubic oversampling), it works as expected:

Then I oversample a.tif with this command:
gdal_translate \
    -r cubic \
    -tr 0.000171662 0.000171662 \
    a.tif \
    b.tif

When openning b.tif in QGIS to generate the hillshade and view it at its native resolution, I expect to get the same result as with a.tif, but here's what I get:

What am I doing wrong? It nearly looks like a nearest resampling is being used (in QGIS I'm using the exact same parameters to generate the hillshades. I'm also getting the same result with a software I am developping).
I also tried in a single step (using -projwin and -tr at the same time), and with gdalwarp, result is always the same.
EDIT: adding some more tests:
If I do this:
gdalwarp \
    -r cubic \
    -te 3.4999998866667914 40.78056666958323 5.968749886666799 43.12434922166657 \
    ../../../emodnet/raw-4326.tif \
    a.tif

gdaldem \
    hillshade \
    -z 1 \
    -s 111120 \
    a.tif \
    a-hillshade.tif

gdalwarp \
    -r cubic \
    -te 3.4999998866667914 40.78056666958323 5.968749886666799 43.12434922166657 \
    -tr 0.000171662 0.000171661 \
    a.tif \
    b.tif

gdaldem \
    hillshade \
    -z 1 \
    -s 111120 \
    b.tif \
    b-hillshade.tif

gdalwarp \
    -r cubic \
    -te 3.4999998866667914 40.78056666958323 5.968749886666799 43.12434922166657 \
    -tr 0.000171662 0.000171661 \
    a-hillshade.tif \
    b-hillshade-2.tif

Then:

a-*.tif have the original resolution, and are displayed correctly in QGIS or Gimp, like in the 1st picture
b-*.tif have the upscaled resolution
b.tif & b-hillshade.tif are not displayed correctly, like in the 2nd picture
b-hillshade-2.tif and the upscaled resolution and is displayed correctly, like in the 1st picture (but I do not need a hillshade tif, I need an upscaled dem)


texts: PostGIS SQL query to find a polygon that intersects with a line

texts:
I need to build a PostGIS query that would select every polygon inside an area the river is flowing through.

I'd better describe with an example below: There is a river (blue line) flowing through several violet polygons (selected in yellow). Both of the yellow polygons are crossed by the green one (labeled Žulová, which is a boundary of a county), but the river actually flows through the only one of them (labeled 2-04-04-044) and the other one is only touched by the end point of the river.

The query I have so far looks like this:

SELECT DISTINCT
        poradi
        ,ST_Distance(ST_StartPoint(ST_LineMerge(toky.geom))
        ,ST_Centroid(pov.geom)) AS dist
FROM
        povodi_iv AS pov
        ,cr_osy_toku_snap_2m AS toky
        ,$table
WHERE
        toky.gid = $gid // river id
    AND
             $table.uid = $id // county id
    AND
             ST_Intersects($table.geom, pov.geom)
    AND
             ST_Intersects(toky.geom, pov.geom)
ORDER BY
        dist ;

But it gives me the result you see on the picture.

Could anyone help me out to build a query that would return only the violet polygons that are both crossed by county boundary and flown through by a river?


texts: What are the other bits in Sentinel's QA60 band?

texts:
According to this page on Google Earth Engine documentation bit 10 is 1 for dense clouds, and bit 11 is 1 for cirrus. The documentation doesn't say anything about all the other 14 bits that can be set, but sends us to this ESA page that explains how QA60 is calculated. On the ESA page, it says that the end result is 0 for clear, 1 for dense clouds and 2 for cirrus. Obviously ESA's output is stored differently than Google's data. From this I take it that Google thought it's a better idea to store QA60 as a bitmask. ESA's output assumes that one pixel cannot be dense cloud and cirrus at the same time, Google's format doesn't.

From my research QA60 is stored on Google Earth Engine on 16 bits. The question is what are the other 14 bits set to, beside 10, and 11, which are dense clouds and cirrus, accordingly. If the other bits are not used, who in their right mind stores 2 bits of data on 16 bits, and starts with bit 10? I could understand 1 byte and start using the first, or the last 2 bits, but why use 16 bits and start in the middle?

Can I safely assume that only 1024 and 2048 values are possible? Or 1536 could also be possible? Or maybe some other values are also possible, which and what do they mean?


texts:
I am able to check if two geometries are intersecting or overlapping through:

 if a.geometry().intersects(b.geometry()):

Is there any method available that tell me if one geometry completely overlaps with another one.


texts: Manually changing elevation data of DEM using QGIS

texts:
I'm looking to do personal modifications to a DEM (Aster global elevation model).
I want to select an rectangular area & reduce its recorded elevation by 20m.
There are other parts on my map I want to raise in elevation, but I think if I learn how to reduce, I will be able to do the process in reverse and increase the elevation.
Is this possible?


texts: Are there any known problems using GPS connection in Android version of QGIS?

texts: I install QGIS both my netbook (windows win7 starter) and My tablet (android 4.0). I can used GPS option in windows version and I can save location points from GPS. However, I could not manage same operation in android verion of QGIS. Is there any problem on android version?


texts: Using internal or external GeoTiff overviews in GeoServer?

texts:
I'm serving some large GeoTIFs using GeoServer. I added internal tilling and overviews. With gdaladdo I have two options:

Adding overviews inside the GeoTIF or;
Saving the overviews in an external file (.ovr), using the read only flag (-ro).

gdalinfo recognizes the external ovr file.

Does geoserver recognizes the external overviews in the .ovr file, or should I insert the overviews in the TIF?


texts: Ordinary Least Squares Checks Prior to Geographically Weighted Regression and what to do when the Jarque–Bera Test is Significant

texts:
I have read that the following 6 criteria need to be met for the OLS model with respect to variables prior to performing a GWR on them: (Is this correct?) and ...

1) Spatial Autocorrelation (my understanding is that spatial auto correlation is acceptable for GWR, since the whole reason for using GWR is to look at local patterns

2) Checking coefficient signs and ensuring that there is a negative for expected inverse relationships

3) Checking redundancy between variables using the VIF (should be less than <7.5)

4) Ensuring variables are statistically significant (p<0.05) and removing ones that are not prior to GWR

5) Ensuring that the residuals of over/under predictions are the result of random noise using the Jarque–Bera (JBT) test, this test result should not be significant

6) Looking at the adjusted R2, a good model is > 0.5

My JBT test is significant, what does this mean in terms of GWR? Can I still proceed or should I transform my data? Use a different model?


texts:
I am trying to plot a Google heatmap based on the JSON objects returned. Here is my JavaScript:

function getHeatMap() {
 heatLayer = new HeatmapLayer({
config: {
    "useLocalMaximum": true,
    "radius": 40,
    "gradient": {
        0.45: "rgb(000,000,255)",
        0.55: "rgb(000,255,255)",
        0.65: "rgb(000,255,000)",
        0.95: "rgb(255,255,000)",
        1.00: "rgb(255,000,000)"
    }
},
"map": map,
"domNodeId": "heatLayer",
"opacity": 0.85
});

$.ajax({
url: "index.aspx/getBusCommuter",
type: "POST",
data: "",
contentType: "application/json; charset=utf-8",
dataType: "json",
success: function (data) {
    var parsed = JSON.parse(data.d);
    $.each(parsed, function (i, jsondata) {
        var coordXicon = jsondata.BusStopX;
        var coordYicon = jsondata.BusStopY;
        var commuterAmt = jsondata.CommuterAmt;
        var data = [
   {
       attributes: {},
       geometry: {
           spatialReference: { wkid: 102100 },
           type: "point",
           x: coordXicon,
           y: coordYicon,
           weight: commuterAmt
       }
   }
        ];

        heatLayer.setData(data);

    });
    map.addLayer(heatLayer);
},
error: function (request, state, errors) {
}
});
}

Basically I tried to set the heat weight by commuterAmt so that it is not just simply plot a heat at each coordinates that I declared. However, it does not work. For every point on the map, it just simply plot a heat without taking account into the commuterAmt.

Attached image shows the heat plotted on each bus stop. Although each bus stop has different commuterAmt but the heat seems all the same.

I wonder is there any alternate way to set the weight as the heatmap used in Google Maps.


texts:
I am trying to extend Symbol class in ARGIS Flex API. this is my widget named TestWidget http://pastebin.com/04w6RSur and here is the extended symbol class http://pastebin.com/bqhDzdSe named CustomSymbol placed in Test package in flexviewer 2.

I had drawn a circle and a line in CustomSymbols sprite like this:

sprite.graphics.beginFill( 0x0000FF,1 );
sprite.graphics.drawCircle(0,0,100);
sprite.graphics.endFill();
sprite.graphics.lineStyle(5,0xFF0000,1);
sprite.graphics.lineTo(mapPoint.x+100,mapPoint.y+100);

when I zoom into map it shrinks its size and when I zoom out it grows accordingly. What I want is make it retain its size whether I am zooming in or zooming out.

Thanks


texts: How to get an area Latlong geometry directly from google map?

texts:
I am new to field of GIS.

I need to know that do we have an option of getting latlong of a building that we see in google earth or any Point of Interest areas.

Something like i create & drag over an area or road(on google map or google earth) and its respective geometrical shape(with Latlong information) gets saved in database or gets popped up.

Is there any tool which does that?

I believe this is what GIS Data Vendors use to create their data.

How can i get geometries of a buildings in a locality boundary or district boundary in few minutes,is it possible?


texts: SVG from QGIS to Illustrator CC gets wrong stroke widths

texts:
When I export SVG files from QGIS using the Print composer the stroke widths become a lot wider when importing into Adobe Illustrator. It seems like the stroke width in the SVG translates differently in Illustrator compared ti Inkscape.

Example

A line in QGIS with the width 0.26 mm is shown as

stroke-width="3.07087"

in the SVG file. And when I open the file in Inkscape, it's still 0.26 mm wide but when I open it in Illustrator it shows the width of stroke as 3.07087 pt.

Does anyone have an idea why Illustrator doesn't "translate" 3.0787 to the right width?

And is it a way of fixing this? I guess I could search and replace in the SVG xml file but I'd rather not..

Does it have to do with the SVG length specs?

  "If a unit is not provided, the length value represents a distance in the current user coordinate system."

Edit: I have now noticed that the mm width is 0.08466669263 times the value in the SVG XML file.

Edit 2: I have noticed now that if I alter or remove the "viewbox" bit in the SVG file, Inkscape also display the line width wrong, so it seems like Illustrator doesn't understand the viewbox and therefore translates the widths wrong!?


texts: GeoDjango - DatabaseError: no such table: spatial_ref_sys with Spatialite

texts:
Im going through the GeoDjango Tutorial. I have loaded the WorldBorders shapefile into Spatialite. I can do spatial queries such contains, intersects, but I have a DatabaseError: no such table: spatial_ref_sys when I'm trying to use GeoQuerySet method such this one:

qs = WorldBorder.objecs.filter(name='Canada')
qs.area()

I have created the database using Spatialite 4.0 GUI. so the SELECT InitSpatialMetaData(); was performed on creation. SELECT * from spatial_ref_sys return all spatial_ref_sys objects.

Can the problem be that the WorldBorders table is in WGS84 and GeoQuerySet method such area() can only be performed with projected SRID?

Thanks for your help

Edit:
The tutorial is the following: https://docs.djangoproject.com/en/dev/ref/contrib/gis/tutorial/

django version: 1.4.3

sqlite3 version is: 3.7.14.1

spatialite version is: 4.0.0

Performing: SELECT * FROM spatial_ref_sys LIMIT 3

-1|NONE|-1|Undefined - Cartesian||Undefined
0|NONE|0|Undefined - Geographic Long/Lat||Undefined
2000|epsg|2000|Anguilla 1957 / British West Indies Grid|+proj=tmerc +lat_0=0 +lo
n_0=-62 +k=0.9995000000000001 +x_0=400000 +y_0=0 +ellps=clrk80 +units=m +no_defs
|PROJCS["Anguilla 1957 / British West Indies Grid",GEOGCS["Anguilla 1957",DATUM[
"Anguilla_1957",SPHEROID["Clarke 1880 (RGS)",6378249.145,293.465,AUTHORITY["EPSG
","7012"]],AUTHORITY["EPSG","6600"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901
"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG",
"4600"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PAR
AMETER["central_meridian",-62],PARAMETER["scale_factor",0.9995],PARAMETER["false
_easting",400000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG",
"9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","2000"]]


texts: Why is FeatureClassToFeatureClass tool Failing on Long cell values (Arc 10.2)?

texts:
I have written a script to convert a CSV into a Feature Class.  The actual line that does the conversion is the very last line in the script in which I call the FeatureClassToFeatureClass tool.

Here is my Code Snippet:

import arcpy
Source_CSV = r"\\int.ec\Data\Stations.csv"
Target_FC = r"Q:\Data\Stations.gdb\Stations"
Description = arcpy.Describe(Target_FC)
Target_FC_name = Description.name
Target_GDB = os.path.split(Description.catalogPath) [0]
Temp_Layer = "Temp_Layer"
SR = arcpy.Describe(Target_FC).spatialReference
arcpy.MakeXYEventLayer_management(Source_CSV, "Lng", "Lat", Temp_Layer, SR, "")
arcpy.FeatureClassToFeatureClass_conversion(Temp_Layer, Target_GDB, Target_FC_name, "", "", "")

I am getting the following error:

File "Q:\scripts\CSVtoASMADPointFeatureClass.py", line 10, in <module>
    arcpy.FeatureClassToFeatureClass_conversion(Temp_Layer, Target_GDB, Target_FC_name, "", "" , "")
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\conversion.py", line 1675, in FeatureClassToFeatureClass
    raise e

ExecuteError: ERROR 001156: Failed on input OID 1263820748, could not write value 'LS_ blue with white trim_ has overflow in Ruisseau Isabelle. Located between "Au P'tit Ruisseau" bar and Bas-Caraquet welcome sign_ just past Rue de Paul. A By-pass was reported at this LS on July 7_ 2012 and a temporary closed status zone was put in place.      ' to output field Comment
    Failed to execute (FeatureClassToFeatureClass).

The value that FeatureClassToFeatureClass could not write is from the 14th row in the CSV. I assume FeatureClassToFeatureClass proceeds sequentially through the CSV, which if correct tells me it had no problem with any of the cell values in the previous lines.

Here are the first 14 rows from my CSV's "Comment" field:

Comment
Foo Bar - overflow to vegetated surface
overflow to vegetated landscape
overflow to field
overflow to field

Flows into adjacent brook.  station is close to lagoon      this lift station serves a few homes - generally will pump 30 mins over two days.  have 2 pumps
Lift station #6 (Rue de l'Ile). No standby power_ dual pumps_ total capacity 500 gal. per min. Overflow pipe to Caraquet Harbour. In case of failure_ alarm system sends signal to pager of LS operator.
LS situated on Chemin St Simon_ No. #7. Located in a small baby barn building off road.
Propane backup; next to lagoon;
October 2015: Back up will be installed at this site this year.   Small grey shingled shack with RV dump station beside LS. Overflow_ near Ruisseau Isabelle. At small picnic park
   October 2015 - lite and buzzer  small station - handles 3 houses and small campground.      Lift station on Rue François Gionet. UTM estimated. This LS pumps at LS Parc Fondateur.      less than 1 %
october 2015 - no lites just pager;  small flow station and no backup   White lift station with siding. Situated at 2264 rue Industrielle. No red warning light_ nor backup power observed.
Blue & white trim lift station on "rue du Lac". New sewage system being installed at time of inspection. LS services 24 houses on Rue du Lac. LS does not have an apparent overflow pipe nor backup power.
LS_ blue with white trim_ has overflow in Ruisseau Isabelle. Located between "Au P'tit Ruisseau" bar and Bas-Caraquet welcome sign_ just past Rue de Paul. A By-pass was reported at this LS on July 7_ 2012 and a temporary closed status zone was put in place.

It has been suggested that the problem is with the single quotation mark in "P'tit", as shown in the error message.

The only problem is P'tit is part of a legitimate place name entered by someone else. I need every character in the comments field to be transferred.

Further, it has been suggested that the problem is not so much with FeatureClassToFeatureClass as it is the CSV converter not escaping single apostrophes.

If this is the problem what coding will force the CSV converter to escape single apostrophes?


texts: Prevent marker from moving in QGIS when zooming

texts:
I was wondering how I can prevent the movement of my markers when zooming out in QGIS.
I have this sample image of Auckland on a larger scale where everything seems to be ok:

Then I zoom out a little bit and everything starts moving:

Also when switching to the print layout everything gets into a new order. Probably I'm just missing something really obvious here...


texts:
I'm trying to display tiles rendered through Mapnik in Mapproxy.

I got it to work correctly for the first zoom level:

But two zoom levels in, the tiles are already combined in a wrong way:

My tiles are served through Mapnik. I followed the following tutorial:

https://www.linuxbabe.com/linux-server/openstreetmap-tile-server-ubuntu-16-04

I can display the tiles correctly through OpenLayers on a slippymap.

Here's my current mapnik configuration:

services:
  demo:
  tms:
    use_grid_names: true
    # origin for /tiles service
    origin: 'sw'
  kml:
      use_grid_names: true
  wmts:
  wms:
    md:
      title: MapProxy WMS Proxy
      abstract: MapProxy Configuration for mapnik

layers:
  - name: mapnik
    title: TMS layer from tiles
    sources: [osm_cache]

caches:
  osm_cache:
    grids: [webmercator]
    sources: [osm_tiles]

sources:
   osm_tiles:
     grid: webmercator
     type: tile
     url: http://10.77.70.19/osm_tiles/%(z)s/%(x)s/%(y)s.png

grids:
    webmercator:
        base: GLOBAL_WEBMERCATOR
        origin: 'sw'

globals:

Some layers further in, everything seems to be okay again:

Is it possible for Mapproxy to show the correct map at all? And if it's not, what are my alternatives?


texts: Charting Normalized Difference of Image Collection?

texts:
I would like to get the normalized difference of each image within an image collection (In this case Sentinel-2) for a specific polygon/geometry and then chart the values calculated for each image.

Is there a way I can do this?

The bands for which I would like to get the normalized difference are B8 and B4.

Sources

Code for the current chart

var filterDate = new Date();
  filterDate = filterDate.getFullYear()
    + '-' + (filterDate.getMonth() + 1)
    + '-' + filterDate.getDate();

  var kcAccum = imageCollection3.filterBounds(Lixin).filterDate('2019-01-01',filterDate);
  var kcChart = ui.Chart.image.series(kcAccum, geom, ee.Reducer.mean());
  print(kcChart);

Results of script

  EDIT

I changed the code to use the map function to compute the normalized difference for each image in the original collection. However, now when I try to chart the computed collection, I get an error message.
The computed image collection is correct. It's just when I try to chart it that it gives problems.

Code Edited

var kcAccum = imageCollection3.filterBounds(Lixin).filterDate('2019-01-01',filterDate).map(function(kcimage) {
    return kcimage.normalizedDifference(['B8', 'B4']).rename('kc');
  });
print(kcAccum);
var kcChart = ui.Chart.image.series(kcAccum, geom);
print(kcChart);

  Error Message


texts: Can I import a kmz file directly from ArcMap 10 to Garmin GPS?

texts:
I currently import a kmz file created in ArcMap 9 (via a very convoluted process using Google Earth) to my Garmin Oregon GPS. I do this so that I can view an aerial photo with town tax parcel map info and other specific features on my GPS. Due to a recent upgrade to Google Earth I can no longer use it on my (very old) computer. We plan to upgrade to ArcMap 10 in the not-too-distant future, and I'm wondering if this will solve my problem, or if I need to purchase a new computer.  Does anyone know whether I'll be able to do this with ArcMap 10 and the Garmin Oregon GPS? Pointers to instructions would be much appreciated!

More about why I do this, and how: I work for a land trust, and often need to find boundaries on individual parcels of land. I sometimes can find either a shape file for the parcels or find the parcel information on Google Earth, but not often, so I go into ArcMap and create a .jpg, then go into Google Earth and add the .jpg as an image overlay. I then find the coordinates of the corners of the .jpg in ArcMap, add these to Google Earth, then click on Save Place As to create a .kmz file which I can move over to my Garmin Oregon device.  It's cumbersome but it works (or did, until Google Earth stopped working on my computer!).


texts:
I have an enterprise geodatabase with about 20 users.  Each user has a corresponding schema. When in Arc Catalog, if only the owned schema box is checked in SQL Management studio I receive a permissions error in Catalog.  However, if both the owned/default schema are selected in addition to the SDE schema, then I have no problem, importing/previewing data.

Additionally, even if SDE is chosen, it grays out, but later becomes un-checked somehow.

For database permissions I have selected create procedure, table, view, and connect.
Roles for each user are configured in SQL Management Studio as db_datawriter.

I am using ArcGIS 10.2.1/SQL Server 2008 R2


texts:
I'm doing a simple query on a points feature class that has a table joined to it. I'm querying for values equal to 0 and my results are incorrect.

No coded domains, this is on ArcMap 10.5. Any explanation for this? Only caught it because I happened to check, but this would be terrible for any analysis.


texts:
Got some problem with Mercator in Geoserver, I tried to display Flight Information Region that has been stored in my Postgre.

The result is more likely drawing using Cartesian projection:

I tried redrawing the FIR based on the same data using WorldWind and got the result:

Here is the configuration of the Coordinate Reference System (EPSG:4326):

And I tried with EPSG:900913 (the result same like the first picture):

Also I tried to change the Native SRS to 900913 but got error:

<ServiceExceptionReport version="1.3.0" xsi:schemaLocation="http://www.opengis.net/ogc http://127.0.0.1:8080/geoserver/schemas/wms/1.3.0/exceptions_1_3_0.xsd"><ServiceException code="internalError">
      Rendering process failed
java.io.IOException
null
ERROR: Operation on two geometries with different SRIDs
</ServiceException></ServiceExceptionReport>

Is there anything wrong with the configuration?


texts: OSM tileserver on Ubuntu 16.04 - renderd doesn't start on reboot

texts:
I have been able to setup a simple OSM tile server on Ubuntu 16.04. I followed this brilliant tutorial and also made some notes of my own. This works alright but whenever I reboot the system renderd doesn't start.

I have been trying to get this working with upstart but that doesn't work on my Ubuntu 16.04 setup.

Here renderd is running fine, I reboot and it's dead. Service renderd start and it's running again.

I have renderd.init in my SystemV init.d directory. Obviously some config must be wrong but I wouldn't know where to start.


texts: Algorithm addgeometryattributes not found

texts:
I am trying to find area and perimeters of Polygon using Pyqgis3. I tried manually first Vector->Geometry Tools->Add Geometry Attributes. By doing so, I get  the area and perimeter column in the attribute table. I would like to automate this process by using python -console inside qgis. I tried the following code, but i get " Algorithm  addgeometryattributes" not found.

Inputparameters={'INPUT':'C:/Demand/war.shp','CALC_METHOD':0,'OUTPUT':'C:/Demand/wardge.csv'}

processing.run('qgis:addgeometryattributes',Inputparameters)


texts: Feature paging with ArcGIS map service vs feature service

texts: I'm about to implement an integration that involves reading all features from a layer on a municipal ArcGIS server.  There are at least two ways to go about this, including use of a Feature Service and a Map Service.  Are the solutions essentially equivalent?  In other words, can I use the Map Service to query and page through results in the same way, or in a way that is equally efficient and problem free, as I can using the Feature Service?  Both solutions appear to involve pulling IDs down and using offsets, there aren't for example 'next' and 'previous' links to utilize, from what I can see.  I am reading documentation on the Map Service presently and see ways forward, but want to see if anyone who has done this can answer this question and possibly provide other guidance as to whether I want to instead go with the Feature Service.


texts: Running LAStools on QGIS 2.8: "missing dependency. lastools folder is not configured"

texts:
To be honest I'm not sure how to run LAStools from within QGIS 2.8.
I'm trying to learn by experimentation.  QGIS shows all the LAStools in the toolbox

If I double click on a LAStool in the toolbox I get the following message:

  Missing dependency. This algorithm cannot run :-(
  LAStools folder is not configured. Please configure it before running LAStools algorithms.

I can't find any LAStools *.exe files anywhere within my QGIS Wien file structure.  That makes me wonder how I could have LAStools show in the toolbox.  Could it be a *.dll somewhere?

What am I doing wrong?


texts:
Using PyQGIS I'm trying to find the chainage point at the start/end of a line based on a chainage points dataset.

I can't figure out how to write the intersect query that buffers the start_point and captures the intersecting chainage points and for now just prints the value.

roadlist=[]
layer = iface.activeLayer()
chainage_layer=QgsVectorLayer(r'Z:\My Drive\Mangoesmapping\Spatial Projects\2019\DSC\126_Road_Updates\Working\Output\DSC_Roads_10m_Chainage.tab', "Chainage", "ogr")
for feature in layer.getFeatures():
    if feature['RoadMntnc'] not in roadlist:
        roadlist.append(feature['RoadMntnc'])

for road in roadlist:
    print (road)
    layer.selectByExpression('"RoadMntnc"= \'{0}\''.format(road), QgsVectorLayer.SetSelection)
    selection=layer.selectedFeatures()
    for feature in selection:
            chainage_layer.selectByExpression('"RoadMntnc"= \'{0}\''.format(road), QgsVectorLayer.SetSelection)
            chainage_selection=layer.selectedFeatures()
            if start_point(buffer(feature.geometry,5,10)).intersects(chainage_selection.geometry()): #based on https://gis.stackexchange.com/questions/256569/select-features-from-another-layer-based-on-a-selection-in-pyqgis/256646
                print (feature2['distance'])
        #aggregate(layer:='DSC_Roads_10m_Chainage.tab', aggregate:='max', expression:="distance" * (attribute(@parent, 'RoadMntnc') = attribute($currentfeature, 'RoadMntnc')),filter:= within(start_point(geometry(@parent)),buffer($geometry,1,10))) # process in model builder -field calculator.

The query needs to run only on the points with the same road name as the line being evaluated.

I get start_point/buffer not defined error.

=====
Here's the updated code

roadlist=[]
layer = iface.activeLayer()
chainage_layer=QgsVectorLayer(r'Z:\My Drive\Mangoesmapping\Spatial Projects\2019\DSC\126_Road_Updates\Working\Output\DSC_Roads_10m_Chainage.tab', "Chainage", "ogr")

for feature in layer.getFeatures():
    if feature['RoadMntnc'] not in roadlist:
        roadlist.append(feature['RoadMntnc'])

for road in roadlist:
    print (road)
    layer.selectByExpression('"RoadMntnc"= \'{0}\''.format(road), QgsVectorLayer.SetSelection)
    selection=layer.selectedFeatures()
    for feature in selection:
        geom = feature.geometry().asPolyline()
        pnt = QgsGeometry.fromPointXY(QgsPointXY(geom[0]))
        f.setGeometry(pnt)
        fbuff=f.geometry().buffer(5,4)
        chainage_layer.selectByExpression('"RoadMntnc"= \'{0}\''.format(road), QgsVectorLayer.SetSelection)
        chainage_selection=chainage_layer.selectedFeatures()
        if fbuff.intersects(chainage_selection.geometry()):
            print (chainage_selection['distance'])

I get ...

  TypeError: QgsGeometry.intersects(): arguments did not match any overloaded call:
    overload 1: argument 1 has unexpected type 'QgsFeature'
    overload 2: argument 1 has unexpected type 'QgsFeature'

Ideally it's best to use an envelope that's the same as the buffer and only select the points within it that are part of the selection_chainage selection to ensure it's the same road. Maybe using something like areaOfInterest = QgsRectangle(fbuff) request = QgsFeatureRequest().setFilterRect(areaOfInterest).setFlags(QgsFeatureRequest.ExactIntersect)


texts: Getting current layer in control event of Leaflet?

texts:
Can I get the current layer in control's event?

    if (obj.overlay) {
        input = document.createElement('input');
        input.type = 'checkbox';
        input.defaultChecked = checked;
    } else {
        input = this._createRadioElement('leaflet', checked);
    }
L.DomEvent.on(input, 'click', this._onInputClick, this);

_onInputClick:function(e){
    var id =e. currentTarget;

}

When user click the checkbox, it will trigger _onInputClick event. This control will be created under a leaflet layer. I want to get the current layer. For example, I can get the current control like e.currentTarget. Is it possible to get the current layer by Leaflet library?


texts:
I am trying to replace values in a single column using an else if statement, but I keep getting an Syntax error in line 2. I'm new to Python so it's probably an easy solution, but I haven't been able to figure it out. I am using ArcMap 10.1.

    def Replace(SrcLevel):
      if SrcLevel == 0:
        return "4-Other"
      elif SrcLevel == 1:
        retun "1-Primary"
      elif SrcLevel == 2:
        return "2-Middle"

    SrcLevel = Replace(!SrcLevel!)


texts: Posting UserCredentials with AsyncTask of ArcGIS for Android?

texts:
My problem is that I want to get JSON data from ArcGIS for Server

For example: http://xxxx.com/arcgis/rest/services/Map/MapServer/legend?f=pjson

I can get and read JSON data from this URL but if ArcGIS for Server has a password I couldn't so , must use UserCredentials but every sample example use UserCredentials for identity task or querytask but I read data from URL.

how can I post UserCredentials with my AsyncTask ?


texts: Adding Layer in current session using ARCPY

texts:
I wanted to add a layer through the toolbox using arcpy.  The below script works fine through the python console in ArcMap, but does not when executed through the toolbox. The layer flashes for a second in the layer management window and disappears. This looks very unlikely and I am not able to figure out how it could not work in the toolbox. Let me know if I had missed something.

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
sourceLayer = arcpy.mapping.Layer('C:/WIP/Shape files/KARTHIK.TEST_TABLE.lyr')
arcpy.mapping.AddLayer(df,sourceLayer,"TOP")
arcpy.RefreshActiveView()


texts: More than one dataset per year in ArcGIS Time Slider?

texts:
Using a Derived/Master Mosaic dataset as a container for source raster mosaics, I have successfully enabled the Time Slider feature in ArcGIS 10.1.  The imagery on display was acquired most years between 2004-2010.   However, in the case of the year 2006 I have more than one dataset; one with leaf-off acquired in early Spring, and one with leaf-on acquired in late summer.

I would like to enable the Time Slider for both raster sets along the same time continuum with the other individual years.  So, in a perfect world I would be able to have: 2004, 2006 (spring), 2006 (summer), 2007, 2008, 2010 on my Time Slider bar.  (I intentionally left out the years for which I have no rasters: 2005, 2009).

However, in order to enable both of the 2006 sets I think I'll have to introduce a Month or even a full date rather than just incrementing by Year.  I'd prefer to not do this because it will cause the user to click 12 times more per year along the slider.

Has anyone come up with a way to enable the Time Slider to include more than one raster per year without introducing month increments?


texts: OpenLayers 3/GeoJSON - Feature 'source'?

texts:
I'm fresh with OpenLayers/GeoJSON...

Anyway, vector layer with features (lines, polygons, points...) and it pulls those data from GeoJSON ("type":"FeatureCollection") file.

Is there a way to see 'feature source' etc. when I click on line, can I find where is the code behind it that 'draws' that line?

My goal is to be able to tell "ok, I have clicked this and it is controlled by this peace of code... or nearby".


texts:
I made a QGIS plugin but I didn't built it with the current structure because I need to keep the dialog window open during execution. So I didn't use a QDialogButtonBox for "Ok" and "Cancel" button, I just used two QPushButton in order to keep the window opened during execution. My problem is when I click on Cancel the dialog closed but the script keep running.

My code is as following :

def __init__(self, iface):
        self.dlg.ANNULER.clicked.connect(self.quitterProgramme)
        self.dlg.OK.clicked.connect(self.exeProgramme)

def exeProgramme(self):
        if programmeExecute == False :
            self.monProgramme()
        else :
            self.quitterProgramme()

def quitterProgramme(self):
        self.dlg.close()


texts: Openlayers 3 - ol.layer.Image rendering on drag'n'drop

texts:
I'm using the lastest version of Openlayers 3, aka 3.15.1.
I want to use

layer = new ol.layer.Image({
    map: olMap,
    source: new ol.source.ImageVector({
        source: new ol.source.Vector(),
        style: layerStyle
    })
});

because i have some layer and a lots a features with hover events, so i wanted to put some layer into an ImageLayer to have better perf.

But when i drag'n'drop my ImageLayer, the part that are out of screen are not rendering when dragging, only when the drop occurs.

Is there something to get a render even when dragging (and animating, interacting, etc...)

i'm ok if i need to modify the openlayers build


texts: Showing menu based on right click on point geometry in JavaScript

texts:
I have implemented a vector layer(using JavaScript) in OpenLayers which has a point geometry.

Now my requirement is to show menu just beside the point whenever we do right-click on the point.

Can we do this in OpenLayers?


texts:
I've been following along with this Ansible Playbook to get PostGIS Tiger data setup:

https://github.com/enigma-io/ansible-tiger-geocoder-playbook

Everything worked great, and it said everything was setup correct with no error messages. The only thing is that it took about 1 hour, instead of the 20-24 hours it said it normally should take (only about 5GB is being used on my volume instead of the 110GB or so it normally does)

So now i'm a bit unsure how to move on - I assume it didn't download all the data I needed (I even checked the config and all the states are set to download).

is there a command I can run or way to proceed to grab all of the tiger geodata/census tracts for each state?  5GB is definitely not enough, right?

Any ideas/help/suggestions?


texts: How to connect PostgreSQL from osm2pgrouting

texts:
I am using OSGeo Live. It is an Ubuntu system. I installed it in my PC. I was following one pgRouting workshop.

https://workshop.pgrouting.org/index.html

In one step, I needed to connect osm2pgrouting with postgis to import data. But, the commands are not working. It is showing following connection fail error. I did not change anything. All the passwords are as default.

$ osm2pgrouting -file "data/sampledata.osm" \
>                           -conf "/usr/share/osm2pgrouting/mapconfig.xml" \
>                           -dbname routing \
>                           -user postgres \
>                           -clean
host=127.0.0.1 user=postgres dbname=routing port=5432
connection failed
$

I tried to fix it. I edited /etc/postgresql/8.4/main/pg_hba.conf, and change the postgres connection method to trust from indent. But, it does not work. The result is same, as before.

# Database administrative login by UNIX sockets
# indent
local   all         postgres                          trust


texts: Retrieving polygon shapes from postgis database and display it to the browser

texts:
I have a codes that will convert polygon shapes that is stored in the database into kml format and count all the rows in it and display it to the browser.

models.py #my model in creating fields in my table

from django.db import models
from django.contrib.gis.db import models as gismodels

class Butuan_Parcel(gismodels.Model):
NewPin = models.CharField(max_length=254)
geometry = gismodels.MultiPolygonField(srid=32651)
objects = gismodels.GeoManager()

def __unicode__(self):
    return self.NewPin;

views.py

from django.shortcuts import render_to_response
from django.contrib.gis.shortcuts import render_to_kml
from butuan_parcel.models import Butuan_Parcel

def all_kml(request) : #functions that convert to kml format
locations = Butuan_Parcel.objects.kml()
return render_to_kml("gis/kml/placemarks.kml", {'places' : locations})

def map_page(request) : #functions that count rows stored in shapefile table
lcount = Butuan_Parcel.objects.all().count()
return render_to_response('butuan_parcel/map.html', { 'location_count' : lcount})

urls.py

from django.conf.urls import *
from butuan_parcel.views import *
from django.views.static import *
from django.conf import settings
from django.contrib import admin

admin.autodiscover()
urlpatterns = patterns('',
(r'^admin/', include(admin.site.urls)),
(r'^kml/', all_kml), #calling the all_kml functions from views.py
(r'^$', map_page)) #calling map_page functions from views.py

map. hmtl

<html>
<body>
var map, base_layer, kml;
function init() {
 map = new OpenLayers.Map ('map');
 base_layer = new OpenLayers.Layer.WMS ("Openlayers WMS", "http://labs.metacarta.com/wms/vmapo", { layers: 'basic'});
 kml = new Openlayers.Layer.GML("KML", "/kml/", { format: Openlayers.Format.KML }); map.addLayers([base_layer, kml]); map.zoomToMaxExtent(); }
 </script>

There should be a total of {{ location_count }} Locations. To see this data in Google Earth, OPen as KML

</body>
</html>

and this is the output in localhost:8000

There should be a total of 105 Locations. To see this data in Google Earth, OPen as KML

now, i want to query the polygon shapes and display it to the browser, not to count anymore.

Is there anyone who knows how to retrieve the polygon and display it to the browser?


texts: How to convert multipatch file in .dwg format to polygon in .shp format?

texts:
I have a multipatch file in .dwg format, which I wish to convert to .shp file format.
When I right click on it to export it to .shp format, then all the data does  not gets converted.

I am using ArcGIS 10.2.2.

I used 3D Analyst -> Conversion -> Multipatch to Footprint, but no result is shown.

How can I convert .dwg multipatch file to .shp file in ArcGIS 10.2.2?


texts: QGIS Welcome page: What is 'pin to list' for?

texts:
In QGIS 3 Welcome page when I right click a project a tiny context menu with 'pin to list' rsp. 'unpin from list' as first entry appears.

'pin to list' produces a little pin near the project name and 'unpin from list' removes the pin:

I consulted the documentation but I could not figure out the deeper sense of these pins.

Does somebody know what these pins are for? Perhaps some kind of favourite projects?


texts: How to performing selection by expression using QGIS

texts:
I am using QGIS 3.8.
I have a grid (hexagon) and a set of points.
I am trying to select the hexagons containing at least 1 point. So far nothing complicated.  Except that sometimes my points overlap 2 hexagons.
So I try to select only the hexagons containing at least 1/2 of the area of the points.
Despite my many attempts, I can't come up with an expression that works.


texts: Does QGIS have WMS-T (WMS with time) support?

texts:
Does anyone know if there's support for WMS-T in QGIS?

Through the app I can add WMS layers, but I can't work with WMS-T (there's no place to set the time parameter).

I been going through the the Python plugins and although there are some that do work with time (Time Manager, Multi temporal and Multivariate data visualisation) they do not support the WMS protocol (they work directly with databases).


texts: Arcgis javascript API & highcharts libraries

texts:
I'm working with jsviewer and JSAPI V:2.8.I try to display a query result in a chart.
what I want to do is :first, I want to  query on my table and  once I got the result liste I have a button that displays the results in a chart.I want to use highcharts librairies.
If there is a code sample or any indication that could help me do that , I will be grateful.


texts:
I'd need to navigate through and mark non-geographical maps.

I used Leaflet.js example to do that (http://leafletjs.com/examples/crs-simple/crs-simple.html), and it's working, but my images would be very large and then have to be served as tiles.

So I installed a local OpenStreetMap server with Mapnik and OpenGIS etc, but I can't seem to find a way to create a .pbf file from my image, hence I can't create nor import the "map" using osm2pgsql...

Is there any way to achieve that ?


texts: What are alternatives to ArcGIS for Desktop and Google Earth for these requirements?

texts:
I currently have a python script that will take in data and output it as a KMZ file, which I then open in Google Earth on my desktop.

The data has two views; a 3D polygon tower (clustered data) and in individual semi transparent circles for a heatmap like view. Additionally, the data can be played over time so that the 3D towers grow/shrink and circles appear/disappear.

The powers that be have told me that the visualization does not "pop" enough for them, and they would like me to look into some free alternatives.

I am afraid that there are no tools besides ArcGIS (which they do not want to buy) that will do what they want, but I figured I would ask anyways. Are there any free tools that meet these requirements?

Ordered by importance:

Can display data over time.
Can read kmz or kml.
Can display 3D polygons and overlay pictures on the map.
Free (or at least cheaper than ArcGIS).

Here are some notes I took while looking at a few free tools:

QGIS (could not get a map to load, could import kml as points)
Open Jump (I think we need our own map data, did not install)
NASA World Wind (Cool, but couldn't get kml to import. Not sure if it has time controls.)
Marble (My favorite, but no time control. Couldn't get kml to input. http://marble.kde.org/kml-guide.php)
MinervaGIS (Have to build from source. Did not install.)
Viking (Worked, no time controls, would not open kmz)


texts: Postgis - Spatial relationships in metres

texts:
I'm trying to perform some spatial relationship analysis using Postgis functions like st_dwithin and st_dfullywithin. I searched and found that for any SRID, using the geography type gives the results in metres which works well for st_dwithin. However, the st_dfullywithin the arguments need to be of type geometry so this approach does not work.

My base SRID is 4283 and I've tried transforming the SRID of the geometry objects to SRID 3577 which still does not work.

As a test, I tried to measure the distance between two points using different SRID and the geometry and geography types.
I believe the units of 3577 is in metres, but when I run the test to get the distance between two points using SRID 3577 it returns the results in degrees, but changing the type to geography returns metres.

I thought changing the SRID to one with metres should return metres? If not, how can you work in metres?

EDIT: The query to get the distance which gives the result in degrees is,

select st_distance(st_geomfromtext('POINT(150.9319747 -33.8482439)', 3557), st_geomfromtext('POINT(151.9319747 -34.8482439)', 3557))

If I use the same SRID, I have to use the following query to get the result in metres

select st_distance(st_geomfromtext('POINT(150.9319747 -33.8482439)', 3557)::geography, st_geomfromtext('POINT(151.9319747 -34.8482439)', 3557)::geography)

My ultimate aim is to fine geoms that are within a given radius of a point,

select
*
from table_1 a
where ST_DWithin(st_transform(a.geom, 3557), st_transform(st_geomfromtext('POINT(150.9319747 -33.8482439)', 4283),3557), 50000);

Inspite of using the different SRIDs, the relationship still considers degrees and not meters.


texts: Featureclass search failing in C++. But working in C#

texts:
I need to perform search on a FeatureClass using ISpatialFilter for a given Point (x,y).
I wrote a program in C# to do that and its retrieving the features as expected. But the similar program in C++ is failing for some feature classes and working for some. I am wondering where it went wrong.

The below c++ code works fine for featureclass cbsa. But doesnt work for BlockGroups.

    //Create Point
    IPointPtr ipPoint = NULL;
    ipPoint.CreateInstance(CLSID_Point);
    ipPoint->PutCoords( dX, xY);

    //Create Envelope
    IEnvelopePtr envelope(CLSID_Envelope);
    hr = ipPoint->get_Envelope(&envelope);

    //Open Feature Class
    IWorkspaceFactoryPtr pWorkspacefactory(CLSID_FileGDBWorkspaceFactory);
    IFeatureWorkspacePtr featureWorkspace;
    IWorkspacePtr workspace;
    hr = pWorkspacefactory->OpenFromFile(CComBSTR(gdbpath,0,&workspace);
    featureWorkspace = workspace;
    IFeatureClassPtr featureClass;
    hr = featureWorkspace->OpenFeatureClass(CComBSTR("CensusBGroups"),&featureClass);

    // Set ISpatialFilter Properties
    BSTR shapeFieldName;
    hr = featureClass->get_ShapeFieldName(&shapeFieldName);

    ISpatialFilterPtr spatialFilter(CLSID_SpatialFilter);
    spatialFilter->putref_Geometry(envelope);
    spatialFilter->put_SpatialRel(esriSpatialRelIntersects);
    spatialFilter->put_GeometryField(shapeFieldName);

    // Search
    IFeatureCursorPtr featureCursor;
    featureClass->Search(spatialFilter, false, &featureCursor); // featurecursor doesnt contain any featueres
    IFeaturePtr feature;

    // Get First Feature
    hr = featureCursor->NextFeature(&feature); //feature is null
    IFieldsPtr totalfields;
    hr = feature->get_Fields(&totalfields);


texts: Unable to download images in Semi-Automatic Classification Plugin (SCP)

texts:
I am using version 3.10, and I am trying to download images directly from the QGIS.
The Semi-Automatic Classification Plugin (SCP) seems to work fine, I can find all the available images.
However, when it comes to download, and I click "Run" and choose the folder, absolutely nothing happens.
Doesn't seem to download anything, but also doesn't give any error.
However, it does some bip sounds like it is trying to connect and doesn't reach it.
Did you have any similar problems? I already reinstalled the Plugin and QGIS.
I followed these settings for MAC, but to no avail.
SCP Plugin for QGIS 3.10 on Mac


texts:
Is there any way I can open Sentinel 2 JPEG2000 file (.jp2) using R? Or more specifically with the {raster} package in R?
I keep getting the Error message that .jp2 is:

  not recognised as a supported file format.

Do you know if this file format will be added to the raster package any time soon, or how to open it in R without converting it to .tif first?


texts:
I have a set of LineStrings that I used to merge with ST_Collect. I want my lines to be curved, so I tried to use St_LineToCurve but my features were still MULTILIGNESTRING. I tried to smooth them on geoserver with stroke-linejoin: round but that didn't work.

I read it's possible to force lines to be curved with ST_ForceCurve. My geometries are now a MULTICURVE type, but GeoServer doesn't seem to handle it :

  Rendering process failed
  java.io.IOExceptionorg.postgresql.util.PSQLException: ERREUR:
  lwgeom_simplify: unsupported geometry type: MultiCurve ERREUR:
  lwgeom_simplify: unsupported geometry type: MultiCurve

Is there a way to curve my multilinestrings?


texts: Trying to maintain spatial reference on newly created numpy array

texts:
I am trying to copy the values from raster into newly created array and then trying to save the array with same spatial reference as raster file

inRas = arcpy.GetParameterAsText(0)

double_layer = arcpy.GetParameterAsText(1)
num_layer = int(double_layer)

output_path = arcpy.GetParameterAsText(2)

# information on input raster
dsc = arcpy.Describe(inRas)
sr = dsc.SpatialReference
ext = dsc.Extent
ll = arcpy.Point(ext.XMin,ext.YMin)

arr = arcpy.RasterToNumPyArray(inRas, nodata_to_value=0)

shape = arr.shape
arcpy.AddMessage("image array created")

freq = numpy.empty(shape=(1, shape[1], shape[2]))
freq_numpy = numpy.uint8(freq)
arcpy.AddMessage("frequency array created")
#count = 1

# check for the 46 layers
# count is giving me total number of pixels

for x in range(0, shape[1]):
    for y in range(0, shape[2]):
        freq = 0
        for z in range(0, num_layer):
            # z is number of layer
            if arr[z][x][y] == 1:
                freq = freq + 1
            freq_numpy[0][x][y] = freq

arcpy.AddMessage("code finish file stored at"+ output_path)

newRaster = arcpy.NumPyArrayToRaster(freq_numpy,ll,inRas.meanCellWdith,inRas.meanCellHeight)
arcpy.DefineProjection_management(newRaster, sr)

newRaster.save(output_path)

The error I get:

ERROR:

Traceback (most recent call last):
  File "D:\PythonScripts\FrequencyBackup.py", line 66, in <module>
    cellHeight = inRas.meanCellHeight
AttributeError: 'unicode' object has no attribute 'meanCellHeight'

Failed to execute (Frequency).
Failed at Thu Dec 17 16:20:53 2015 (Elapsed Time: 4.75 seconds)

When I try to access inRas.meanCellHeight in Python console of ArcGIS, it returns a value. How can I access it from within my script?


texts: How to select only polygons that overlap with another layer?

texts:
I have two shapefile layers created from different thresholds from a predictive model (Pink=80; Yellow=95).

Is there a way to create a new layer of only those pink polygons that contain yellow polygons and effectively erase the pink polygons that do not contain yellow?

If I use Geoprocessing - Intersect I will end up with just a copy of the yellow layer.


texts:
I want to show cluster of point data on map using OpenLayers. I have been following this example https://openlayers.org/en/latest/examples/earthquake-clusters.html but the problem is that want a look like dot density map like each state has only one dot and dot size varies according to count.
var style_ch = new ol.style.Style({
    pointRadius: "${radius}",
    fillColor: "#FF3300",
    fillOpacity: 0.8,
    strokeColor: "#66FF00",
    strokeWidth: "${width}",
    strokeOpacity: 0.8,
    label: "${getLabel}",
    fontSize: "12px",
    fontFamily: "Courier New, monospace",
    fontWeight: "bold",
    labelOutlineColor: "white",
    labelOutlineWidth: 2.5
}, {
    context: {
        width: function (feature) {
            return (feature.cluster) ? 2 : 1;
        },
        radius: function (feature) {
            var pix = 2;
            if (feature.cluster) {
                pix = Math.min(feature.attributes.count, 7) + 2;
            }
            return pix;
        },
        getLabel: function (feature) {
            if (feature.cluster) {
                if (feature.cluster.length > 1) {
                    return feature.cluster.length;
                }
            }
            return '';
        }
    }
});

layer = new ol.layer.Vector({
    source: new ol.source.Cluster({
        distance: 30,
        source: new ol.source.Vector({
            format: new ol.format.GeoJSON(),
            url: function (extent) {
                return "http://localhost:8080/geoserver/db/ows?" +
                    "service=WFS" +
                    "&version=1.1.0" +
                    "&request=GetFeature" +
                    ..... +
                    "&outputFormat=application/json";

            },
            strategy: ol.loadingstrategy.bbox
        }),
    }),
    // style: style_ch
    styleMap: new ol.style.Style({
        "default": style_ch,
        "select": {
            fillColor: "#8aeeef",
            strokeColor: "#32a8a9"
        }
    })
});

But it is simple cluster closed points and each state has more than one value according to the distance.
I want to show dot per every state on the basis of feature count of the vector layer. Here I got a problem that the styling doesn't apply and still due to distance property map show many points of cluster instead of one. I want my output result looks like this:


texts: I have two layers: polygon and raster. Could you tell me how to count the number of pixels which intersect with polygon layer?  Raster is too large to creating fishnet or convert to vector (time consuming). Also zonal statistic doesn't work because it gives me the average pixel value and not their number. I would ask for help, any program.


texts: Convert geometry type from Polygon to Multipolygon using FME?

texts: I don't want to alter the geometry, just change the geometry type from Polygon to Multipolygon, in order to comply with a constraint.


texts: Limiting cloud estimation of Landsat8_SR over ROI in GEE

texts:
I intend to create an image collection of Landsat-8 surface reflectance not TOA that has 100% free cloud over my Region of Interest (ROI). Not surprisingly the code I have found in the forum ( Filter Landsat images base on cloud cover over a region of interest ) did not work for me.
Based on https://gis.stackexchange.com/users/51905/intotecho, I tried to clip the collection first by var clippedLan8 = lan8.clipToCollection(ROI); where lan8 is my collection and of course ROI is the geometry defining my area, but it does not work.
So I am asking: 1) Is the code mentioned applicable for surface reflectance (SR) or not? If not, is there a way to achieve the goal in GEE? 2) Why doesn't the code for the clip work?


texts: How can I print a bunch of Arcmap mxd's all at once?

texts: I have a folder with a dozen Arcmap compositions (.mxd files) I need to print. They're all ready to go, I don't need to do any editing. Do I really have to open up each one, wait to refresh, and press print?


texts: Connection QGIS 3.4.5 / PostgreSQL 9.3 with role group

texts:
I created 2 connection roles with PostgreSQL. Then, to manage more easily the privileges, I created a group role includind these roles (inherit).

I created a trigger which track data changes. I used the session_user function to know who made this action (Insert/ Update/ Delete) but it got me the group role login. I want to get the user (ie. role)  login and not the group role login.

I read that there is another function named current_session but it not get the user login.

Is there a another function I didn't find ?


texts: Sentinel-2 Query scihub apihub wget missing results

texts:
When I perform a search in Firefox using
https://[username]:[password]@scihub.copernicus.eu:443/dhus/search?q=footprint:"intersects(41.9000, 12.5000)" AND ((platformname:Sentinel-2))

I only get to see 10 of the 443 results with no option to see the others, same in IE.
When I use wget to get results as XML it's the same, I get an XML output with 10 results using:
wget  --continue --config wgetrc "https://[username]:password]@scihub.copernicus.eu:443/dhus/search?q=footprint:\"intersects(41.9000, 12.5000)\" AND ((platformname:Sentinel-2))" -O c:/downloads/test6.xml

How do I get it to produce an XML output with all 443 results?
A manual search in the hub does allow me to see all foubd images


texts: Selecting values with date-format by using SelectLayerByAttribute_management with input Parameter

texts:
I am working on a Tool which selects records depending on the input value.
Works fine, till I want to select date-records. I used

def buildWhereClause(self, table, field, value):

    # Add DBMS-specific field delimiters
    fieldDelimited = ap.AddFieldDelimiters(table, field)

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

    # Add single-quotes for string field values
    if str(fieldType) == 'String':
        value = "'%s'"% value
    elif str(fieldType) == 'Double':
        value = float(value)
    elif str(fieldType) == 'Date':
        value = datetime.strptime(value,'%d.%m.%Y')

  That's my idea of the quote which I added to accept date-format in the inputValue. But when I am typing in for example 07.01.2019 (a date which exists in my table), I get this error: ExecuteError: ERROR 000358: Invalid expression
  Failed to execute (SelectLayerByAttribute).

    # Format WHERE clause
    whereClause = "%s = %s" % (fieldDelimited, value)
    return whereClause

def execute(self, parameters, messages):
    """The source code of the tool."""

    inputFeature = parameters[0].valueAsText
    inputFieldname = parameters[1].valueAsText
    inputValue = parameters[2].valueAsText
    outputFeature = parameters[3].valueAsText

    whereClause = self.buildWhereClause(inputFeature, inputFieldname, inputValue)

    ap.SelectLayerByAttribute_management(inputFeature,"NEW_SELECTION", whereClause)

    return

I am wondering if my understanding of the function of this method is basically wrong, or if it's just an issue with the syntax...


texts: Extract all pixels values from GeoTIFF with Python

texts:
I would like to recover the values of the pixels of a GeoTIFF of a single band. First I get the x and y coordinates and I would like to use them to access the pixel value. Here is my code:

   raster = rasterio.open('.../g.tif')
   br_pixel = (raster.width * raster.height)

   liste3 = []
   i = 0
   while i <= nbr_pixel:
      (x,y) = raster.xy(i,i)
      listX.append(int(x))
      listY.append(int(y))
   i = i + 1
   l = list(zip(listX,listY))
   liste3 = raster.sample(l)
   print(liste3)

   <generator object sample_gen at 0x7f7ef5821c50>


texts:
I have many raster tiles (in tif format, SR WGS84) and, since all together are very heavy, I want to unify them in a geopackage.

I've tried to do that in QGIS, performing the following steps:

converting all the raster in PNG format,
creating a geopackage from the Browser panel (Create Database > No geometry)
simply dragging them into the geopackage created ---> Import succesful

but: the original rasters have a single band, while once added to the geopackage they have 4 bands.

Is there a way to create a Geopackage that contains all my raster tiles (about 850 tiles) maintaining only the 1 original band?


texts: Getting "Null" for the joined columns values after joining a shapefile and CSV file in QGIS

texts:
I am trying to join one CSV file and one shapefile.

This is the shapefile attribute table:

and this is the CSV file.

I am getting this as the result:

Does anyone know why I am getting null values in the TyincomeUsa_In... column?

Edit:
What all options do I check in this for the given problem in the case of this CSV data


texts: I grouped the layers in ArcGIS Desktop and published them to my Portal for ArcGIS, but when I tried to create the Map Application, I see the individual layers. They are not in a group.


texts: How to find the country in which a point is located?

texts:
I have a PostGis database created with osm2pgsql. I am looking for a way to add the country name in which a point is located to the point row. It must work also for countries with enclaves, such as France (Llívia). In this case, osm2pgsql creates several polygons from the border relation: a polygon with the outer border and and the enclave, and for each of the both countries a polygon with the enclave only. A simple ST_Within or ST_Contains does not work in this case. I found I way to do the job:

update planet_osm_point p set country=
       (SELECT distinct poly.name
         from planet_osm_point point inner join planet_osm_polygon poly
         on (st_within(point.way, poly.way)
            and (st_nrings(poly.way)=1 or not st_within(point.way, st_makepolygon(ST_InteriorRingN(poly.way,1))))
            and not exists (select * from planet_osm_polygon other
                              where other.admin_level='2'
                                    and other.name=poly.name
                                    and st_contains(st_makepolygon(st_exteriorring(other.way)), poly.way)
                                    and other.way_area>poly.way_area))
            where  poly.boundary='administrative' and poly.admin_level='2' and point.osm_id=p.osm_id)
      where p.name is not null;

This seems a bit complex to me and takes a lot of processing time. Is there any simpler and faster way to do this?


texts: Creating new features from common part of other features in same layer using QGIS

texts:
I have many overlapping features in one layer and I want to get common parts of it (like in intersection function, when using two different layers). Moreover I want to sum one of features' attributes and assign it to common parts. Like in the above picture.

Is there any way to do it easy and/or quickly?


texts:
So my situation:

I've got an array of 4 microphones, 3 in a plane and one above out of plane. These microphones are connected to a micro controller which accurately time-stamps a sound when it arrives, in microseconds. Since the time the signal was made is unknown, these 4 time-stamps just give 3 time-difference-of-arrival measurements, i.e., how long after the first detection the other 3 microphones detected the same signal.

From these 3 TDoAs I'm trying to pin point the location of the sound source. It is possible for me to add another microphone if needs be but research has indicated that it should be possible using only 4 microphones.

Obviously the positions/coordinates of the microphones are known, as well as the speed of sound.

Is this possible without knowing the time the sound signal was made?

I've spent quite a lot of time researching this and have tried writing codes to perform triangulation and multilateration but to no success.


texts:
I have been trying to solve the problem to get a small portion of a larger GIS map into my AutoCAD Civil 3D for a site plan. I currently obtained a large map from the county as a .gdb file and I have successfully imported it.

I just however want a small portion and when I import into Civil 3D I want it to have elevation references which it currently has in the attribute table. I am currently trying to use GRASS plug in to accomplish this.

The steps I have taken so far have been:

Made a new shapefile layer
performed a toggle edit and made a polygon around the area
used processioning tools to make a clip of the area that I want the elevations to show.

I am left with the portion I want however in GRASS it cannot find the layer (clipped layer) to take the elevations off when I use the 3D command.

Can anyone suggest how to associate the elevations with my vector layer or give me an easier suggestion on how to import the data I need into Civil 3D with the elevations given in the contours?


texts: gdaldem hillshade using command to calculate hillshade of whole folder

texts: I need to calculate hillshade of folder using gdaldem hillshade which has 60 files in it. How can i do it using command prompt window.


texts: Unable to rename attribute table columns using QGIS

texts:
I'm using QGIS and trying to remove the "ACS block group" prefix from each field in this attribute table.  When I try to rename my attribute columns it says "Failed to rename field to ...Is this field unique?"
It is a unique field, so what else might be happening?


texts: Raster multiplication script not executing

texts:
I built this code to multiply rasters but it will not execute, it gives me the following error: ERROR 000732: Input Raster: Dataset in_raster1 does not exist or is not supported

# Import arcpy module
import arcpy
from arcpy import env
from arcpy.sa import *

arcpy.CheckOutExtension("3D")
arcpy.CheckOutExtension("spatial")
arcpy.env.overwriteOutput = True
arcpy.env.workspace = outRas = "P:\del\New Folder\new"

# Local variables:
in_raster1 = arcpy.GetParameterAsText(1)
in_raster2 = arcpy.GetParameterAsText(2)
in_raster2 = arcpy.GetParameterAsText(3)

# Process: Raster Calculator
outRas = Raster("in_raster1") * Raster("in_raster2")* Raster("in_raster3")

Any help?


texts: SEXTANTE Saga-grid normalisation and KBS

texts: I tried to use the saga tool grid normalisation. The result is ok but the KBS has changed. Is this an error?


texts:
Regarding on the previous question asked https://gis.stackexchange.com/a/303372/119333, solution from best answer is implemented and its working. I am rendering tab panel inside popup overlay and you can click on tabs without click bubbling further on map. Now I added button inside my tab panels and when i click on that button (yellow) it bubbles to map and executes map.on('click',function()), ignores if statement despite fact mouseOver variable is true.

Code:
Button rendered inside table row in tab panel:

BROJ UGOVORA	<button id="brugovoratrigger" onClick="pozoviUgovor(this)" style="background:yellow">'.$row['broj_ugovora'].'</button>

JS:

var ug=document.getElementById('brug');

ug.addEventListener('mouseover',function(){
map.getInteractions().forEach(function(interaction){
interaction.setActive(false)
mouseOver=true
   })
})

map.on('click',function(e){
if(mouseOver){ return }
else {
map.removeOverlay(popup)
var coor=e.coordinate
...

Here is picture of what I am talking about:


texts: Is the software used in the "seabed UFO" video GIS software? If so, what is it?

texts:
The icons almost make it look like the software is ArcReader or an ArcEngine based, but it's hard to tell from the video.  I supposed it could not be GIS software at all.

Is anyone familiar with this software?

Here's the video link:
http://www.youtube.com/watch?v=jwHE5816_XY&feature=player_embedded#t=56s


texts:
I have made the following SQL view on Geoserver:

SELECT *
FROM areas
WHERE landuse = '%landuse%'

This works perfectly fine when I specify the landuse in the URL, for example 'grass'. However when I want the query to return all features regardless the type of landuse I do not know what I have to enter in the URL.


texts: QGIS can't connect to PostGIS but pgAdmin4 works

texts:
I have a PostGIS DB at a remote server and I can connect to it from my desktop using pgAdmin4. However, when I try to connect though QGIS (v3.4) I get an error:

FATAL: no pg_hba.conf entry for host "10.1.<complete ip here>..., user "car_r" database "car", SSL off

The funny thing is that my computer is listed on the server pg_hba.conf or else I would not be able to connect using pgAdmin4.

Any idea what could be wrong?

Adding more info

My PostGIS server is setup to accept only secure connections (hostssl in the pg_hba.conf). When I force QGIS to use SSL, (SSL mode: require) I get a "SSL error: certificate verify failed" error.

How can I make QGIS recognize my server certificate?


texts: Triangulation method to create multipatch objects (decomposition of polygon surfaces)

texts: Which triangulation method or algorithm does arcscene use to create multipatch parts (patches) as you see in the attaching picture. The algorith decompose the polygon area into a set of triangles using vertices of the polygon object. But which method is it? It is not delanuay or ear clipping in my opinion but which one I can not figure it out any help??