texts: WMTS in leaflet (from standalone geowebcache)

texts:
I have a standalone geowebcache, I want to use it to serve wmts to a leaflet map but I find it hard to understand what url to use.
Displaying the layers in QGIS work fine when I add the server using the url http://localhost:8080/geowebcache/service/wmts/
But how do I find out what the url should be for leaflet?
If I check getCapabilities i find http://localhost:8080/geowebcache/service/wmts/rest/test_layer/{style}/{TileMatrixSet}/{TileMatrix}/{TileRow}/{TileCol}?format=image/png
But how do I implement that?
I am using default style. I have created a custom TileMatrixSet called "EPSG:3857"
I tried this without luck
L.tileLayer('http://localhost:8080/geowebcache/service/wmts/rest/test_layer/EPSG:3857/{z}/{y}/{x}?format=image/png')

Previously I have successfully displayed layers in leaflet using tms from geoserver using the structure below, but for some reason i get 400 error now.
L.tileLayer('http://localhost:8080/geowebcache/service/tms/1.0.0/test_layer@EPSG:3857@jpeg{z}/{x}/{-y}.jpeg')

Really scratching my head here because I try to understand where I should look and and how I can implement this properly? Is it possibly to use wmts with leaflet? I found a wmts plugin for leaflet https://github.com/alexandre-melard/leaflet.TileLayer.WMTS but is that really the only/recommended way to do it?
EDIT:
This url is working
L.tileLayer('http://localhost:8080/geowebcache/service/wmts?layer=test_layer&style=&tilematrixset=EPSG:3857&Service=WMTS&Request=GetTile&Version=1.0.0&Format=image/jpeg&TileMatrix=EPSG:3857:{z}&TileCol={x}&TileRow={y}')
But I dont understand how I can get there from looking at getCapabilities alone? I found out by looking in the network tab for from geowebcache/demo. But how I know how to from `
http://localhost:8080/geowebcache/service/wmts/rest/test_layer/{style}/{TileMatrixSet}/{TileMatrix}/{TileRow}/{TileCol}?format=image/png

to
http://localhost:8080/geowebcache/service/wmts?layer=test_layer&style=&tilematrixset=EPSG:3857&Service=WMTS&Request=GetTile&Version=1.0.0&Format=image/jpeg&TileMatrix=EPSG:3857:{z}&TileCol={x}&TileRow={y}

that I dont understand?


texts: Changing (flipping) line direction in QGIS

texts:
I have a river network in QGIS. I draw opposite direction of one of the rivers. Now I want to change the direction of that river.
How can I do that in QGIS 3?
(As I know it is flip in ArcGIS)


texts: How do I specify the position, style and units of the ScaleBar in pyqgis?

texts:
I would like to add a scale bar to my composition using the python console and then export it to a raster image.
I am following the MapRendering tutorial with the "output using map composer " approach
however,the Scale Bar item example do not have many features.

I would like

a scale bar with 2 ticks, one at the beginning, one at the end
units in km
eventually round it to the nearest km
place it on the upper right of my figure.
set a transparent background

Just like in the desktop QGIS version.
What are the options ?
So far I have found

item = QgsComposerScaleBar(c)
item.setStyle("Line Ticks Up") # optionally modify the style
item.setComposerMap(composerMap)
c.addItem(item)

Will the scale bar be "at scale" ?
Thank you for your help

I cannnot find anything from that page
ScalBar class


texts: Change value of marker cluster using folium

texts:
I would like to know if it was possible to change the value displayed in the marker clusters wich, by default, display the number of markers it contains.

I think this customization of marker cluster (depicted here) as implemented in leaflet would fit my needs. But I don't know if it is possible to do so with folium.

My goal is to display an average value contained in the children markers of the clusters.


texts: ERROR 000271: Cannot open the input datasets Failed to execute (CompositeBands) (arcpy)

texts:
I am trying to stack bands 3 and 4 together for numerous landsat scenes, and I am trying to do it using a loop and arcpy.CompositeBands_management.  The code I am using is:

import arcpy
#set environment
arcpy.env.workspace=r'F:\Sheyenne\Atmospherically Corrected Landsat\hank_masked\NDVI'
#path to band 3
band3 = arcpy.ListRasters("*band3*")
#path to band 4
band4 = arcpy.ListRasters("*band4*")
#loop that matches raster names and only executes process for those that match
for raster in band3:
  rasterName = raster
  stack = "" #initialize your variable to be populated in the inner loop
  for band in band4:
    if band.startswith(rasterName.split("_")[0]):
      stack = band
      break #break out of the inner loop - we have a match
  if stack != "": #just to make sure we have something
    #generic naming
    outraster = rasterName.replace(".tif", "_stacked.tif")
    arcpy.CompositeBands_management('rasterName;stack',outraster)
print('Done Processing')

but this returns the error:

ExecuteError: ERROR 000271: Cannot open the input datasets
Failed to execute (CompositeBands).

I am thinking it has something to do with this line:

arcpy.CompositeBands_management('rasterName;stack',outraster)

If i take the ' out of this line it still won't work though

Edit:

If it helps at all, when I print out the variables rasterName, stack and outraster they are populated as I would expect


texts:
I would like to convert a shp file to map file. I know should use these steps:

shp -> osm -> pbf -> map

For step 1 I use a python script : ogr2osm.py

For step 2 I use osmosis, but I get an ERROR:

SEVERE: Thread for task 1-rx failed
org.openstreetmap.osmosis.core.OsmosisRuntimeException:
Node -13963 does not have a version attribute as OSM 0.6
are required to have. Is this a 0.5 file?


texts: Creating thematic map using OpenLayers and GeoServer?

texts:
I am still a beginner in learning GeoServer and OpenLayers.
I tried to create thematic maps of the map that I publish in GeoServer. I used OpenLayers to show the map in web browser. When I want to create thematic maps based on data from the attributes that I have, I find the color constraints of thematic maps can not arise.

What to make thematic maps should be used SLD in GeoServer?

Here this code:

<!DOCTYPE html>

<meta charset='utf-8' />
<title>indoprop</title>


texts: How do I use scipy.ndimage.gaussian_filter1d on an arbitrary axis?

texts:
I have a GDAL raster that looks like this:

And I would really like to blur this raster along an arbitrary transect.

I see that scipy.ndimage has a gaussian filter that allows me to blur the raster along a given axis, as long as the axis is valid given the raster's numpy matrix.  Doing so with the default axis of -1 (and sigma=2) looks like this:

This raster is pretty close to what I'm looking for, I'd just really like to have this filtering not strictly happen along the horizontal axis, if that's possible.

So my questions are thus:

Is it possible to filter the matrix along an arbitrary axis, such as an axis that is not necessarily along the vertical or horizontal?
Is it possible to do so just with Numpy/Scipy and GDAL?


texts: Loading a GeoJSON in QGIS using Python but the layer can't be edited

texts:
def loadLayer(self):
    layer = self.dlg.listWidget.currentItem().text()
    url = self.getRoute(layer)
    data_req = requests.get(url)
    data_layer = data_req.json()
    with open('xx.geojson','w+') as f:
        json.dump(data_layer,f)
    layer = QgsVectorLayer('xx.geojson',layer,' ogr')
    QgsProject.instance().addMapLayers([layer])

I'm getting data from a server in the form of a feature collection.  I write it into a GeoJSON file then load it in QGIS. The layer loads properly but if I try to edit some features it is not possible.


texts: PostGIS raster connection is not closed by MapServer

texts:
I'm using PostGIS raster in MapServer like this:

TYPE RASTER
DATA "PG:port=6432 dbname='mydb' user='postgres' password='1234' schema='public' table='truemarble' mode='2'"
PROCESSING "CLOSE_CONNECTION=ALWAYS"

My problem is that the database connection is not closed after a request. The PROCESSING directive should ensure the disconnection and this works fine with PostGIS vectors.

Any ideas?

I've tested with:
python_mapscript in a mod_python application
php_mapscript in a php-fpm application (Fastcgi)

MapServer 6.2 and 6.4, GDAL 1.10, Ubuntu 12.04 and 14.04


texts: Importing representations from Layer File?

texts:
I'm trying to incorporate representations as part of my data model. I have set up representation rule for my feature class, which is part of a feature dataset.

The way I have the data categorized is by unique value field. There are four categories - A(105 records), B (35 Records), C (54 Records) and D(38 Records). I have also renamed the labels. I've created the representation in ArcMap and I can see this fine in the FGDB.

I'd like to import this representation symbology into a clean copy of the same layer, which resides in a different location. I am using 'Add Representation' tool to achieve this, however the output is not correct.

I have saved the representation FC to a layer file, which acts as the 'Import Rule Layer' in the 'Add Representation' toolbox. However the output's representation FC symbology, although classified to 4 categories from the initial representation, aggregates all values into the first symbol - A, which now contains all 232 records and these are symbolised according to category A. The representation does not see any records in classes B,C or D, however I can see these remaining categories in the representation tab.

I'm not sure why this is occurring.


texts: Forcing y-axis range in ui.Chart.image.series()

texts:
I am harmonizing Landsat 8 and Sentinel-2 datasets for agricultural mapping purposes. When plotting the results using ui.Chart.image.series function, one graph shows the y-axis starting not at 0. For comparison purposes I need the y-axis on the charts to start at zero. How to adjust the below command to force y-axis to start at zero?
// Plot the fitted model and the original data at the ROI.
print(ui.Chart.image.series(fittedHarmonicboth.select(['fitted','NDVI']), roi, ee.Reducer.mean(), 10)
    .setOptions({
      title: 'Harmonic model: both S2 and L8 values',
      lineWidth: 1,
      pointSize: 3,
}));


texts:
I am trying to read a csv file straight into a set of new feature classes within a new file geodatabase, but I am getting the following stack trace when it fails:

  Attempted to read or write protected
  memory. This is often an indication
  that other memory is corrupt.
  ESRI.ArcGIS.Geodatabase at
  ESRI.ArcGIS.Geodatabase.IFeatureDataset.CreateFeatureClass(String
  Name, IFields Fields, UID CLSID, UID
  EXTCLSID, esriFeatureType FeatureType,
  String ShapeFieldName, String
  ConfigKeyword) at
  WaterEditorExtension.Model.SurveyUtilities.StreamWriting(String
  pointtype, String outfile, String[]
  record, IWorkspace workspace) in
  C:\Documents and Settings\lbadgerow\My
  Documents\Visual Studio
  2008\Projects\WaterEditorExtension\WaterEditorExtension\Model\SurveyUtilities.cs:line
  206.

I believe that it is related to a schema lock, but I am able to create a feature dataset which is confusing.  (In proper stacktrace form, below is my workspace factory call, and my where I'm building the FC)

        public IWorkspace CreateFileGdbWorkspace(String path, String jobnumber)
        {
            Type factoryType = Type.GetTypeFromProgID("esriDataSourcesGDB.FileGDBWorkspaceFactory");
            IWorkspaceFactory workspaceFactory = (IWorkspaceFactory)Activator.CreateInstance(factoryType);

            IWorkspaceName workspaceName = workspaceFactory.Create(path, jobnumber, null, 0);

            IName name = (IName)workspaceName;
            IWorkspace workspace = (IWorkspace)name.Open();
            return workspace;
        }
...
IFeatureWorkspace featwork = (IFeatureWorkspace)workspace;
featwork.CreateFeatureClass("test", outfields, CLSID, null, esriFeatureType.esriFTSimple, "SHAPE", "");

EDIT 17 March 2011:  has anyone attempted to Thread a component (extension) of ArcMap?  I'm thinking that I could release the schema lock that is preventing this from moving on to the FC creation phase of this project, but I can't seem to get ArcMap to remain stable when I create a new thread to build the gdb.


texts: No Grass plugin in plugin manager in Qgis 2.14 on ubuntu 16.04

texts: I have installed QGIS 2.14.15-Essen and GRASS GIS 7.2.1, but  there is no option to add GRASS plugin in Plugin Manager


texts: Copying features from 'in_memory' to 'in_memory' using ArcPy?

texts:
I came across this post and found out about Copy_management().

Is it possible to copy an in_memory feature to another in_memory feature?

I tried to test it using the following code, which fails:

shp = r"D:\test\somefile.shp"

bufferOne = arcpy.Buffer_analysis(shp,r"in_memory\bufferOne",'1 meters')

copyOfFirst = arcpy.Copy_management(bufferOne,r'in_memory\BufferOneCopy')

Trying the exact same thing with absolute paths, instead of in_memory works fine. Therefore I assume you cannot copy from in_memory to in_memory, or am I missing something? Would there be a way to achieve the same effect?


texts: Where to get raster files of Earth, united states

texts:
Possible Duplicate:
  Open access repository of general GIS spatial data?

I'm brand new to GIS programming so please excuse my lack of knowledge.

I am developing a GIS web app using MapGuideOS to host raster images (base maps?) and PostGres for some vector data. I've been reading a lot of documentation and books but some things are still a little confusing...

Where can I get some free base maps for the united states (Does that mean I need shapefiles)? I want to host something similar to the basic view of google maps (state boundaries, county boundaries and roads), but only for the United States. I know I can use google maps and draw my vector data on it, but I would like to explore options for hosting it myself, before I commit to GoogleMaps/Bing/Yahoo.


texts: ESRI JavaScript HeatmapRenderer class and Print widget

texts: I’m working on an application using ESRI’s HeatMapRenderer class from the JavaScript API.  The issue that I am running into is that the Print widget fails if the heat map renderer layer is on the map. If I do not have the heat map layer on the map the Print widget works as expected.  I was wondering if anyone else has run into this issue, and if so what sort of solution was found.


texts: Is there any reason that case of characters in Proj4 strings would affect their interpretation? Before parsing a Proj4 string, can I lower case all the characters first?


texts: Extend polygons to match another polygon in ArcGIS

texts:
I am struggling with this task :(
I have various layers, all having a common outer shape (see the blue polygon in the picture).

On of those polygons however, has a different shape.

I somehow need those polygons (the colorful ones) to extend, so they match the outline of my other polygon. The values do not really matter much, as long they are plausible relating to the neighboring cells.

I tried with a buffer, but it didn't really work. Manual reshaping is also no option, as I have hundreds of those areas.

Any idea?


texts: Diverse values to one field for soil analysis

texts:
I have soil analysis data in an Excel table. For one analysis point there are values at three depths. Besides the Excel file I also have a shapefile with the points to which I want to add the soil data. What is the common way to add three values in an field of one row or similar?

Example (Excel table):

ID  point  depth   pH-value   n-value  n-value   n-value   ~40-more-values
1   1      0-20    4          .....
2   1      20-40   5
3   1      40-60   7
4   2      0-25    4
5   2      25-50   5....

This could look like the Excel file or like the following structure in a shapefile:

ID  point  depth              pH-value   n-value   ....
1   1      0-20;20-40;40-60   4;5;7      ...
2   2      ...

I think this is not nice for doing queries. But when I make one column for each depth and soil value the table would have about 120 columns.

Also when having three rows for one sample point and each depth (like in the Excel table), then I actually have three points at the same coordinates in my shapefile. Is this a good idea when I put the shapefile into a WebGIS app and want the user to get the information about the soil-sample-point by a mouse click?

Is there a common way to represent this kind of data?


texts: Removing small spaces (slivers) between polygons?

texts:
I have a set of polygons representing administrative districts.  The districts were digitized by hand, and there are small spaces (slivers) between polygons where the polygons should be touching.
How can I have the polygons snap together, removing the spaces?


texts:
I want to define a featureId fiter for wfs layer like below:

wfs = new OpenLayers.Layer.Vector("WFS Vectore", {
        strategies: [new OpenLayers.Strategy.BBOX(), saveStrategy],
        projection: new OpenLayers.Projection("EPSG:4326"),
        protocol: new OpenLayers.Protocol.WFS({
            version: "1.1.0",
            srsName: "EPSG:4326",
            url: "http://localhost:8080/geoserver/iran/wms?service=WFS",
            featureType: "population",
            featureNS: "http://iran.kadaster.org",
            geometryName: "the_geom"
        }),
        filter:
            new OpenLayers.Filter.FeatureId({
                fids: ['population.913', 'population.912']
                //type: ?????
            })

    });

I don't know to what set 'type' variable in filtering option?


texts: What SRS does leaflet use when it creates new markers?

texts:
I'm using leaflet 0.7. I'm wondering what spatial reference system is used when a new marker is created? Also, how does this impact my use of maps? (I'm using leaflet .7 because of dependency requirements from other plugins). I'm using ESRI base map tiles from the esri-leaflet plugin, which I think are EPSG:3857. I have data in a postgis database in EPSG:4269. When I use the latlong from a new marker in leaflet, I get errors when trying to use that to compare it to spatial data in my database - IE find nearest geometries. I keep reading on the internet that leaflet uses 3857. However, when I transform the point from 3857 to 4269, my spatial queries don't work. However, when I do something like this, and assume the point is in 4326 in leaflet and then transform to 4269, it works - ST_Transform(ST_SetSRID(ST_MakePoint(-96.78427, 32.78299),4326),4269)).

So -

Is this the case that a new marker in leaflet has coordinates stored in EPSG:4326?
If I add a map tile layer like ESRI that is EPSG:3857, or some other coordinate system, even though I transform the coordinate clicks from the markers for use in the DB, does this mean that the points don't quite line up as they should on the map?


texts: How to slice a list based on row values?

texts:
I'm working with python 2.7 (quite a beginner).
I have a list, and I have a feature class with :

A = a field that contains values that I would like to use as the start of the slicing

B= a field that contains values indicating the number of items I need to take into my slicing.

Actually, I would like to do something like this :

list = [1,2,3,4,5,6,7,8,9,10,11,12]

rows = ap.SearchCursor(fc, fields = "A; B")

for row in rows:

    a = row.getValue("A")

    b = row.getValue ("B")

    newlist= list[a: (a+b)]

But its not working, I got an error "slice indices must be integers or None or have an index method"


texts:
I am new on Python and using Python within ArcGIS and a bit stuck hoping someone could point me in the right direction.

I want to remove the 6 last characters if the field is longer than 6 characters.
I'm using the script in the field calculator.

def doThis(fld):
    val = 0
    if len(fld) > 6:
        val = fld[len(fld)-6]
    return val
    else:
        val = fld
    return val

Dothis( !Asset_ref!)

With the code above I get a SyntaxError: Invalid syntax (line6)


texts: How to draw (or add) a line indicating a specific latitude/longitude?

texts:
I'm trying to add a single line that show a specific latitude  (e.g. the Tropic of Cancer) I don't want a gird but only a single line even in a vector layer.How can I do that I've try with pyQGIS console but seems I can't...

I would like to measure angle with that line so I'd like to have it in a layer.
However I found useful a method to use it in print composer, but my main goal is to have it as line in QGIS.

For the moment I'm using Google Satellite open layers with a EPSG:3857.


texts: Problem transforming projection on outlying features from shapefile

texts:
I am fairly new to GIS data and as such please feel free to ask for extra data I might not have provided.

I have some regional data from Natural Earth, more specifically this: http://www.naturalearthdata.com/downloads/10m-cultural-vectors/10m-admin-1-states-provinces/ (States and Provinces, first download).

This uses a longlat-projection, which I wish to turn into a miller projection. However, several of the regions get "mangled" by the transformation. It seems to only affect those sitting at the border of the map (far west / far east), more specifically:

(In the format "Country - Region name")
Antarctica - ATA-00 (Antarctica aggregation)
Brazil - Goiás
Fiji - Northern
Russia - Chukchi Autonomous Okrug
Zambia - Eastern

If I manually go in and edit the nodes touching the edge, it seems to fix the problem. However, I am trying to automate the process, so I can automatically generate new data from the shapefile when a new release is made.

Is there any smart way of overcoming this issue? The error presents itself as if the feature is trying to combine all the polygons into one, drawing ugly lines across the entire map to connect them. An example can be seen here: http://i.imgur.com/bzJ5V.png (taken from QGIS). The result is the same whether through a script (using GDAL) or through the QGIS program.


texts:
Looking for a free spell check for ArcGIS 10.1.
Someone wrote a VBA for pre-10 here.
Is there a way to implement it in 10.1, perhaps writing it in python instead of VBA?


texts: Leaflet - get bounds of all layers, not just vector

texts: I need to get the bounds of all layers in a Leaflet map, not just vector layers. For vectors, I add them to a L.featureGroup, which has the method getBounds(). How do I include raster layers as well?


texts: Clipping raster with different UTM zones using rasterio?

texts:
I have 2 Landsat images, with an overlapping area:

The one on the left is in WGS84 / UTM32N (EPSG:32632), the one on the right is in WGS84 / UTM33N (EPSG: 32633).
I have a polygon in WGS84 / UTM33N (EPSG: 32633). The polygon is located in the overlapping area and I would like to clip each one of the rasters according to my polygon.

Theoretically I should have a problem with the left raster because it is not in the same CRS as the raster, but using clipping tool in QGIS 3 I was able to clip both rasters and the results is satisfying. Each new clipped raster got the same CRS of its 'parent' raster.

The problem rises when I try to perform the same task using rasterio:

src = r"c:\rasters"
fields_geojson_path = r"C:\Ran\tries\fields\tryField_UTM33N.geojson"
with fiona.open(fields_geojson_path, "r") as shapefile:
            geoms = [feature["geometry"] for feature in shapefile]
            geoms = list(geoms) #this is because the original file is geojson, if it was shp this needs to be deleted

for i in os.listdir(src):
    if i.endswith('.gz'):
        pass
    else:
        file_list = []

        for j in os.listdir(os.path.join(src,i)):
            if 'stack' not in j:
                if j.endswith('.tif'):
                    date = j.split("_")[3][2:]
                    new_tif = os.path.join(src,i) + "\\"+ date + "_LS8_" + j.split("_")[8].split(".")[0]+".tif"

                    with rasterio.open(os.path.join(src,i) + "\\" + j) as src3:
                        out_image, out_transform = mask(src3, geoms, crop=True, all_touched=True)
                        out_meta = src3.meta.copy()

                    out_meta.update({"driver": "GTiff",
                                     "height": out_image.shape[1],
                                     "width": out_image.shape[2],
                                     "transform": out_transform})

                    with rasterio.open(new_tif, "w", **out_meta) as dest:
                            dest.write(out_image)
src.close
dest.close
print "finish"

The clip with the UTM33N raster works fine, but when the iteration get to the UTM32N raster I get:

ValueError                                Traceback (most recent call last)
<ipython-input-61-bf98345b5021> in <module>()
     17
     18                     with rasterio.open(os.path.join(src,i) + "\\" + j) as src3:
---> 19                         out_image, out_transform = mask(src3, geoms, crop=True, all_touched=True)
     20                         out_meta = src3.meta.copy()
     21

C:\ProgramData\Anaconda2\lib\site-packages\rasterio\tools\mask.pyc in mask(*args, **kwargs)
     10     import warnings
     11     warnings.warn("Deprecated; Use rasterio.mask instead", DeprecationWarning)
---> 12     return rasterio.mask.mask(*args, **kwargs)

C:\ProgramData\Anaconda2\lib\site-packages\rasterio\mask.pyc in mask(raster, shapes, nodata, crop, all_touched, invert)
     66     if rasterio.coords.disjoint_bounds(source_bounds, mask_bounds):
     67         if crop:
---> 68             raise ValueError("Input shapes do not overlap raster.")
     69         else:
     70             warnings.warn("GeoJSON outside bounds of existing output " +

ValueError: Input shapes do not overlap raster.

How can I do that with Rasterio?

Can I force it to perform it somehow?

My developing platform is Jupyter notebook, Python 2.7 on Windows, so I would prefer some guidance in Python libraries, preferably Rasterio.


texts:
I have to trasform a script (which I didn't write myself) from QGIS to ArcMap.

In the script a filter by expression is used. The function has the input inl (= the layer with the data that should be filtered) and filter (a string e.g. 'HQ_SZEN<=30' / '"HQ_SZEN"<=30'). This filter is a user input and in QGIS one has to use the second form. I think ArcMap requires the same input but it didn't seem to work so I tried to use arcpy.AddFieldDelimiters with the first form like mentioned in example 4 here https://pro.arcgis.com/en/pro-app/arcpy/data-access/searchcursor-class.htm. To use this I had to split the attributename from the actual expression and my reasoning was that there would always be a strange character between the name and the expression. So that is the reason why the arcpy code is kind of ugly with the chars_allowed_in_layername.

I was thinking of solving it by using arcpy.da.SearchCursor with a where clause but I can't get the filter to work. As filter i use 'HQ_SZEN <=30', the output I get is:

HQ_SZEN
<=30
"HQ_SZEN"<=30
No filter - filtering skipped

Main code:

fifeatures = addAttrFactory().FilterInput(fl, filter)
        try:
            arcpy.AddMessage("Filter returned {0} rows".format(len(fifeatures)))
        except:
            arcpy.AddMessage("No filter - filtering skipped")

QGIS code:

    def FilterInput(self,inl, filter):
        ftext = ''
        if not filter or not filter.strip():
            #print "FilterInput()..."
            featureIterator = inl.getFeatures()
        else:
            ftext = filter
            #print "FilterInput(%s) ..." % (ftext)
            featureIterator = inl.getFeatures(QgsFeatureRequest().setFilterExpression ( ftext))

        #print "FilterInput(%s) done" % (ftext)
        return featureIterator

My ArcPy code:

    def FilterInput(self, inl, filter):
        ftext = ''
        if not filter or not filter.strip():
            featureIterator = arcpy.da.SearchCursor(inl, '*')
        else:
            chars_allowed_in_layername = "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ_0123456789"
            for char in filter:
                if char not in chars_allowed_in_layername:
                    end_layername = filter.index(char)-1
            name = filter[:end_layername].strip(' ')
            arcpy.AddMessage(name) #with input HQ_SZEN<=30: HQ_SZEN
            filter_expr = filter[end_layername:]
            arcpy.AddMessage(filter_expr) #with input HQ_SZEN<=30: <=30
            ftext = u'{0}'.format(arcpy.AddFieldDelimiters(inl, name)) + filter_expr
            arcpy.AddMessage(ftext) #with input HQ_SZEN<=30: "HQ_SZEN"<=30
            featureIterator = arcpy.da.SearchCursor(inl, '*', where_clause=ftext)

        return featureIterator


texts: Error ET Geowizard Clean Polygon() takes no arguments (3 given)

texts:
I want to clean polygon feature classes with ET Geowizards (clean polygon) because I cannot find any tool within ArcMap (ArcInfo) that does the same. I have really small gaps within the polygons (but no overlaps), which should be removed.
The process is running well (the most time) in ArcMap. I am using a large python script with different steps and it has worked as well the first time but since a few days the script fails at the ET Tool with the CMD message "ET Geowizard Clean Polygon() takes no arguments (3 given)". The tool does not start the process, it fails at the beginning.

I did not change the input data, code or software! The ET logfile does not contain any information, only "Error". I am using the same data input for testing the whole time, is has just a few MB. I tried the tool via CMD on different machines, in ArcMap 10.1 and ArcMap 10.2. Always the same.

I have installed the last update of ET Geowizards (11.2 for ArcGIS 10.2, 14 Dec 2014)
Python version 2.7.5 64 bit (not running with 32 bit because of memory (2 GB memory is NOT reached during the ET process)
ArcMap 10.2.2

The data input is a feature class within a FGDB, output the same, Fuzzy tolerance 0.5 meters. The code was running well... but now?

(myWorkspace = arcpy.GetParameterAsText(1) (FGDB), LN_STREETS_DISSOLVE is an output feature class that is created before)

try:
    CleanPolygonResult = myWorkspace + os.sep + "final_streets"

    arcpy.AddMessage("Clean Polygons")
    arcpy.ET_GPCleanPolygon(LN_STREETS_DISSOLVE, CleanPolygonResult, "0.5")
    arcpy.AddMessage("Clean Polygons done")

except Exception as ex:
# If an error occurred
    arcpy.AddError("Clean Polygons failed: " + str(ex))
    sys.exit()

Any ideas why the Tool fails? I would also be glad if I could solve the problem without the ET Geowizards, only using ArcMap tools but it has to run within the python script, not by manual editing :) Nevertheless I am interested in the problem of the ET Tool...


texts: Is there a way to assign a symbology automatically based on the file name in QGIS?

texts:
In QGIS is it possible that when I create a shapefile and give it a file name e.g. Site, that I can assign symbology to that automatically?

I have a set symbol for Site and many other features but would love it if I could auto assign the symbology based on the file name.

I know this is possible in ArcMap but am trying to find out whether this can be done in QGIS.


texts: Working with QGIS Zanzibar. I have a data set where many values are NULL (or blank), I wish to filter the data to show only those with NULL. In the query builder I use "FieldName" = NULL, but this removes ALL the data... Why. The data I'm trying to filter on is a date.


texts: Is it possible to remove small polygons from Census Block shapefile?

texts: I just pulled up a shapefile that plots census blocks (TIGER product) using Quantum and I noticed that there are well over 1000 small polygons that aren't census blocks that are (I assume) used to buffer between census block polygons and small slivers of open area.  Is there any way I can remove these irrelevant polygons so that I am left with a clean "grid" of census blocks that are all adjoining without these small polygons?


texts: Converting Excel into shapefile (keeping encoding and data type)

texts:
I am trying to create a "perfect" shapefile from the Excel data. "Perfect" means that it should match several factors: proper encoding and original data type.

Excel file basically contains information about point objects and includes five attributes with required types.

Attribute     |      Type       |      Description
--------------------------------------------------------------------------------
FFA_NR        |      integer    | Object serial number
FFA_TXT       |      text       | Object serial number but with leading zeros
Name          |      text       | Object name
X             |      real       | X-coordinate in Gauß-Krüger CS
Y             |      real       | Y-coordinate in Gauß-Krüger CS

where an X,Y-tuple is based on DHDN / 3-degree Gauß-Krüger zone 5 EPSG:31469

Data sample looks as following (demonstrated in the image below).

My workflow (similar to this Importing Spreadsheets or CSV files.) so far is

In Excel, Data > Save As...
I am trying with .xlsx, .txt(MS-DOS) or .csv(MS-DOS)
In QGIS everything visualized nicely in the sense that it is placed correctly. I have done it with several approaches, Layer > Add Layer > Add spreadsheet layer and Add Delimited Text Data
Afterwards, in QGIS, I am saving the opened point layer as a shapefile with Save vector layer as...
Then in the Attribute table, I do have several issues with data encoding and data type

1. Data Type

When I import Excel data into QGIS I am partially losing my initial data types. For instance, instead of the desired 035, I do have 35. Obviously, it is about the import data type.

I have seen similar questions on this topic.

Changing CSV layer attribute value from text to number in QGIS?
How to keep leading zeros in csv file when working with QGIS

Keeping leading zeroes in QGIS in big files
How to format labels with leading zeroes?
Joining .csv file to a layer in QGIS [closed]
How to Specify Data Types of CSV Columns for Use in QGIS

However, I can not understand how to create the file with .csvt extension. What is the workflow? How to create it and where it should be stored?

Of course, in Attribute table I can create a new column, appropriated with a text type, and then proceed with a Filed Calculator, using for instance lpad("FAA_Nr", 3, '0').

2. Encoding

I am working with German text, that contains letters like ä, ö, ü, ß. When in QGIS, I apply any of the encoding formats (System, UTF-8, ISO-8859-1, ISO-8859-15), none of my specific characters can be visualized.

Check the following graphics for UTF-8 and ISO-8859-15 accordingly.

I have found related topics,

Exporting from Excel to csv and import in QGIS retaining encoding/special characters with workaround?
Converting Hebrew Excel file into Shapefile in QGIS

However, I did not discover anything spectacular in the context of my issue. I am guessing that the problem is more with the Excel rather then QGIS, is not it?

Or maybe I apply the encode formatting on a wrong step. When should it be done?

Target: optimizing the efforts required for the data input in QGIS.


texts:
I know there isn't a way to export a map from ArcMap to PDF and keep the hyperlinks, but is there a way to add hyperlinks to the map in Adobe Acrobat? I have a simple map of a few states and I am wanting each state to link to a webpage. I know how to get to seeing the attributes in Adobe Acrobat, but can't figure out if I can add a hyperlink there.

I think this is different to Exporting PDF with Hyperlinked Attributes from ArcMap? which instead explains how you would have to go through the model tree to find the link. I am asking how to create a link on the states on the pdf, so that you wouldn't have to search through the pdf to find the links.


texts:
I followed this example of putting a d3 graph in a leaflet popup. The code returns with div.node(). What exactly is a node? I would like the Leaflet Popup To display this custom HTML, and then have the d3 chart below it:

var content = '
'+feature.properties.SubDivName+'
' + document.getElementById("sasi").innerHTML + document.getElementById("chart")

I'd like to keep the chart() function separate from this HTML. But when I try and combine the output from the chart function and this bit of HTML, the d3 chart does not display. How can I get the chart() function to return my D3 graph as just another HTML element I can append to my 'content' variable?


texts: When I copy and paste a feature, why does some of the attributes fields show null?

texts:
I work on some maps representing administrative-territorial units in Roumania during 1990-2016. So, for every year that are associate with modifications in those units, I must realize a map.
I downloaded a shp for 2016, and then I merged the features for representing the situation of 1990.
In 2000, it has been established two new units. So, I copied these two features from 2016 shp and paste on 1990.
But, after I did it, some attribute fields for these two copied features appear as null, even in the first layer (shp-2016), the values are shown.

How to resolve this problem?


texts: Optimize Earth Engine cloud cover area calculation for Sentinel-2

texts:
I am trying to generate cloud free Sentinal 2 images, I have a basic cloud detection code working but it runs very slowly (most likely due to nested map and reduce features). I was wondering if there was a way to do this better?

function ratioRelevent(im){

// Area of the image tile
  var ftImg = ee.Algorithms.GeometryConstructors.Polygon( ee.Geometry( im.get('system:footprint') ).coordinates() )
// Intersection of the image and the ROI
  var inter = roi.intersection( ftImg, ee.ErrorMargin(0.5) )
// Percentage of the image covered by the ROI
  var coveragePer = ee.Number(inter.area()).divide(roi.area())

  var roiImgCloud = im.select(["cloud_mask", "B2"])
        .reduceToVectors({
          reducer: ee.Reducer.mean(),
          geometry: roi,
          scale: 100,
          geometryType: 'polygon',
          eightConnected: false,
          labelProperty: 'constant',
          }).filter( ee.Filter.eq('constant', 0) )

  var addArea = function(feature) {
     return feature.set({areaHa: feature.geometry().area(ee.ErrorMargin(0.5))});
     }

  var roiImgCloudArea = roiImgCloud.map(addArea)

// the sum area of the cloud
  var areaSum = ee.Number(roiImgCloudArea.reduceColumns({
      reducer: ee.Reducer.sum(),
      selectors: ['areaHa']
      }))

  var areaSumNum = ee.Dictionary(areaSum).get('sum')

  var CloudPerc = ee.Number(areaSumNum).divide(inter.area())
//the cloud percentage of an image WRT the ROI
  var CLOUDY_PERCENTAGE_ROI =  ee.Number(areaSumNum).divide( roi.area()) ;

  im=im.set('CLOUDY_PERCENTAGE', CloudPerc)
  im=im.set('cover_ROI', coveragePer)
  im=im.set('CLOUDY_PERCENTAGE_ROI', CLOUDY_PERCENTAGE_ROI)
  return im
}


texts: Georeferencing help (GeoTiff into MapWindow6)

texts:
I'm a newbie with MapWindow6 and I would like to figure out a couple of things.
1 - Whenever I drag and drop a GeoTiff into the map (MapWindow6) I can not see that image at its place, I mean that Image is not georeferenced, so if I put another GeoTiff the latter overlaps the first one.
With MapWindow4 it was quite simple after creating an instance for the grid raster.
2 - Is there any particular way to implement DotSpatial in my project? I'm using Visual Studio 2005.
Any tips would be really appreaciated.
Thanks a lot guys,
Carmelo


texts: Explaining difference between shx and shp files of shapefile?

texts:
I am searching for a more in depth explanation of the differences in shp and shx files to little avail. I mean beyond the 'shp contain geometry - shx contains an index of the geometry'.

The reason I ask, is because while working in QGIS recently, I have made two observations that led to questions in my mind around the exact differences in these file extentions:

QGIS can open and display both shx and shp, and each file appears to be more or less identical in their output (display),
but not exactly so - I have noticed that sometimes the matching shx/shp files display slightly 'off-kilter' relative to each other. It doesn't appear to be a projection issue, they simply don't draw in the exact same location as each other.

These observations made me curious as to why these differences in display exist, and why QGIS can open and operate the shx in the same manner as the shp, when previously my understanding was that the shp is the 'master' file if you will, but requires .dbf and .shx to function correctly as a single, whole entity.


texts: I read the discussion (Openstreet with leaflet not appearing on some desktop and all mobile browser) but I don't know where to find this code and fix it on Joomla 3. I used a plugin (Fields - GMapFP OpenStreet Map) that works fine on desktop browsers but doesn't show nothing on mobile (android/chrome).   Link: my website


texts:
I set up my QGIS plugin development environment in Eclipse with PyDev using the instructions.

I've got the autocomplete working for qgis and qt4 objects:

but I can't figure out how to get the code complete working for my form controls:

Notice how I only get autocomplete for functions I've already used.

I realize this is tricky since the control is defined in the UI file. I just wanted to check if anyone else had solved this.


texts: Which API should we use for printing a PageLayout?

texts:
We are building an standalone ArcGIS Engine application and we are going to add support for printing.

There are several APIs for printing a PageLayout in ArcObjects, in ArcGIS 10 there are at least three different ways to do this:

The new PrintAndExport class, in
the documentation they seem to
promote this API Conceptual help -
Printing maps.
It's also  possible to print a
layout using the PrintPageLayout method PageLayoutControl Sample: Printing with the PageLayoutControl
But also Output method on IActiveView can be used.

Can someone give recommendations on which API we should choose?  Is there any known pitfalls in these APIs?

Update:
Found one pitfall during my research:

When working with ArcGIS Server layers, only PrintAndExport draw the patches/swatches on the legend when printing the layout.


texts: How to get the Nth LineString from a MultiLineString in PostGIS?

texts: I need a function that gets the Nth LineString from a MultiLineString in PostGIS, somthing like ST_PointN that gets the Nth point from a LineString.


texts: Source for North American public boat and kayak launch sites?

texts: Does anyone know of a good source for North American public boat and kayak launch sites for mapping in a GIS. I would prefer a web service.


texts:
In python how do I subprocess.call the gdaldem module for color relief?
In the command line I would type this:

gdaldem color-relief path-to-raster path-to-colorfile path-to-destination-raster -alpha


texts: Plotting via GeoPandas but only one point shows. What's wrong with my code? Map and points are in same projection

texts:
I've been stuck at this for two weeks and have tried already solutions online. I am using Python 3.8.

from shapely.geometry import Point
import geopandas as gpd
from geopandas import GeoDataFrame
import pandas as pd
import os

os.chdir(r'path')
df = pd.read_csv('emscPhilippines2008to2020.csv', delimiter=',', skiprows=0, low_memory=False)
geometry = [Point(xy) for xy in zip(df['Longitude'], df['Latitude'])]
crs = "epsg:32651"
gdf = GeoDataFrame(df, crs=crs, geometry=geometry)
basemap = gpd.read_file('PH_provs.shp')
gdf.plot() # POINTS ONLY
gdf.plot(ax=basemap.plot(figsize=(17,15)), marker='o', color='red', markersize=15);

Both are in EPSG:32651. Even if the CSV is imported in 4326, same result.

Result

Points only


texts: Getting items from OrderedDict to parse addresses?

texts:
My task is to take a full address and parse them out.  I found a module called usaddress ( https://pypi.python.org/pypi/usaddress ) that will take '12345 Old Hillbilly LN Lot 7 Oklahoma City, OK 78945' and output the OrderedDict below.

How do I go about assigning and populating them to the correct field?

(OrderedDict([('AddressNumber', u'12345'), ('StreetName', u'Old Hillbilly'),
('StreetNamePostType', u'LN'), ('OccupancyType', u'Lot'), ('OccupancyIdentifier', u'7'),
('PlaceName', u'Oklahoma City'), ('StateName', u'OK'), ('ZipCode', u'78945"')]), 'Street Address')

This is what I have so far just to get it to populate one field called AddressOut but it errors.

import os
import usaddress
import arcpy
from collections import *

# Get feature class through parameter
fc = arcpy.env.workspace = arcpy.GetParameterAsText(0)

field1 = "FullAddress"
field2 = "AddressOut"

with arcpy.da.UpdateCursor(fc, field1,field2) as cursor:
    for row in cursor:

    row.setValue(field2, row.getValue(list(usaddress.tag(field1))) # field2 will parse field1
    cursor.updateRow(row)

New code:

import usaddress
from collections import OrderedDict

# Get feature class through parameter
fc = arcpy.env.workspace = arcpy.GetParameterAsText(0)
field1 = "FullAddress"
field2 = "StreetName"

rows = arcpy.UpdateCursor(fc)

for row in rows:
    addr_dict = usaddress.tag(row.getValue(field1))
    x,y = addr_dict
    #arcpy.AddMessage(addr_dict)
    a = x[field2]
    arcpy.AddMessage(a)
    row.setValue(field2, a)
    rows.updateRow(row)
del row
del rows


texts:
One thing that I really like about leaflet popups is that they have the autoPan option, which means that if I select a feature programatically in my JavaScript, the map will scroll to fit it, without me having to pass the map around everywhere in my code.

Recently, we changed our spec to get rid of popups from our map. But I still want the autoPan feature! I'm tempted to create an "invisible" popup with code like this, so the popup never shows up, but is open just long enough to trigger the autopan.

  layer.bindPopup('', { closeOnClick: false })
  layer.on('popupopen', () => layer.closePopup())

But this seems really hacky. Does anyone know of a better way to autopan to a feature, without having a reference to the map?


texts: Trilateration algorithm for n amount of points in R?

texts:
I read the page Trilateration algorithm for n amount of points and I have to do similar kind of trilateration in R.

I am not that expert in R, can anyone help, if there is similar function/package like mathematica in R?


texts: Copying Table of Contents from QGIS into Excel?

texts:
Is there a way to copy and paste the contents of the Table of Contents in QGIS into an Excel spreadsheet?

For example, I have many layers (groups of vector shapefiles which are underneath various 'groups' which I would like to maintain/keep these categories when transferring the file names across to the sheet) which I would like to preserve when copying and pasting.

It is not the attribute table itself, and I can't just copy and paste the legend of a print out map as I need to be able to have this as text.


texts: Deleting features with arcpy on ArcGis SDE 10.0 SP5

texts:
I am new to ArcGIS and arcpy and am learning largely through trial and error. Anyway, the overall task is to delete, edit and add rows to feature classes.

I thought it would be good to start with delete - it seemed the simplest. In the desktop I can change to a version, start an editing session and delete the object.
However, when I try the script below, it fails with a database error. On the desktop I don't seem to be able to delete unless I do 'Start Editing', but I can't find that in arcpy for 10.0. I see over the years other posts but nothing quite like this. Is this something that I should be able to do with arcpy?

arcpy.env.workspace = myconnect
versions = arcpy.ListVersions()
versionName = "montgs1"
layer='layer'

arcpy.MakeFeatureLayer_management("ARCFM.DistributionMain", "layer")
arcpy.ChangeVersion_management(layer, 'TRANSACTIONAL', user + "."+versionName,'')

qry ='"OBJECTID" in (3142049,4228649)'

arcpy.SelectLayerByAttribute_management(layer, "NEW_SELECTION", qry)

numsel =int(arcpy.GetCount_management(layer).getOutput(0))
if numsel > 0:
    print str(numsel) + ' features selected'
    arcpy.DeleteFeatures_management(layer)

Incidentally, when I use that criteria in the gui, it selects two items, however in the script numsel's value is 1.

Anyway, when the delete happens I get this:

  (, ExecuteError('ERROR 999999: Error executing function.\nUnable to create logfile system tables. User perhaps lacks permissions or resources to create tables\nUnderlying DBMS error [ORA-01455: converting column overflows integer datatype\n] [ARCFM.DistributionMain]\nFailed to execute (DeleteFeatures).\n',), )

Why does it work on the GUI but not with arcpy?


texts: Is there any way to force a loop to go until the last list object in R?

texts:
I am running a loop in R to rasterize a list of data frames and for comparison purposes I have to keep the same list length in the output.

I am facing some problems with NULL objects on the list. I managed to handle NULL objects inside the list with tryCatch, however when I have null objects at the end of the list (e.g. last 3 objects) they are not included in the output of the loop, so the output list (fry) get shorter than the original list (fr.fr). See the output below.

I would like to know if there is an efficient way to handle this problem in R and to ensure the same list length. I also looking to more straightforward approaches.

Here is my code:

# Raster
st <- raster(xxx)

# List of data frame to Raster
fry <- list()

for (i in 1:length(fry.fr)){
  tryCatch({
    fry[[i]] <- rasterize(fry.fr[[i]], st, fry.fr[[i]]$AREA, fun = sum, na.rm=FALSE)
  }
  ,error= function(e) {cat("Error", "\n")})
}

Here is my output:

  length(fry.fr)
   336

  length(fry)
   333


texts: Find flipped lines in ArcMap geometric network?

texts:
I have been working on a project tracing lines in a geometric network.  I have had to reverse the line (flow) direction on many lines to make the traces work.  Now I am at a point where I have to update the master dataset, and so make those same line reversals in the database.

My thought is that I'd run a small arcpy script to check the to/from XY for every line against the to/from of the matching lines in the master dataset - any that have the to XY matching the master from XY (and vice-versa) would be recorded for flipping.

Although I think this script shouldn't be too difficult to write, I am now wondering whether there is a tool already available that will identify the flipped lines for me?  I haven't been able to find one in ArcToolbox yet, but I may be looking for the wrong thing.

Note there are other geometry changes in my dataset that I don't want to find, just identical matching lines (matched on an ID) where the only difference is direction.  Lines where XY is different I want to ignore.

For example, the follow lines showing line direction

The lines in the blue area I want to take note of, the only difference is the line direction.  The line in the red area I do not want recorded as the line has actually moved.

Is there a tool in ArcGIS Desktop (10.5) that can identify lines that have been flipped, but ignore other geometry changes?  If so, which one?

I have just found the "Detect Feature Changes" tool, which I am testing right now (it is very very slow).  I am unsure whether this is going to find just flipped lines, or everything.  Of course I may have set the parameters incorrectly.

The Detect Feature Changes tool finished processing - took close to 20 minutes to run - but only found 9 features where the line direction had changed.  I am expecting at least 10000 flipped features.


texts: Using R data with ArcPy to make histogram in ArcMap?

texts:
I am writing a tool in python which passes parameters to a R script, which in turn creates a DBF file. I read this DBF and want to make a histogram out of the data inside ArcMap as one of the results of the tool.

The code snippet is as follows:

import arcpy

simloglkhd = "H:\\CIV\\KULLDORFF\\DATA\\PLOTDATA.dbf"
inputTemplate = "C:\\Program Files (x86)\\ArcGIS\\Desktop10.1\\GraphTemplates\\Default.SimplifiedChinese.tee"
output_graph_name = "VerticalBarGraph1"
verticalBarGraph_bmp = "H:\\CIV\\KULLDORFF\\DATA\\plotdata1.bmp"
# Create the graph
graph = arcpy.Graph()

# Add a vertical bar series to the graph
graph.addSeriesHistogram(simloglkhd, "simloglkhd",10)

# Specify the title of the left axis
graph.graphAxis[0] = "Frequency"

# Specify the title of the bottom axis
graph.graphAxis[2] = "log(lambda)"

# Specify the title of the Graph
graph.graphPropsGeneral.title = "Monte Carlo Distribution of Lambda"

# Output a graph, which is created in-memory
arcpy.MakeGraph_management(inputTemplate, graph, output_graph_name)

# Save the graph as an image
arcpy.SaveGraph_management(output_graph_name, verticalBarGraph_bmp, "MAINTAIN_ASPECT_RATIO", "600", "375")

I am having the following difficulties:

I am not able to find a template for this histogram. In the GraphTemplates folder there is PopulationPyramid, but it is still not working. Is there somewhere I can find templates? The create template operation gives a text file instead of a .tee file.
I get an error in the line where I addSeriesHistogram saying "Field '10' is not in the dbf"(something like this) while it is only the countBin. Can someone tell me where I can find good examples of the usage of this function?

My dbf file is simple
ID(1-999) &
simloglkhd(floating point numbers).

I really do not know how use this MakeGraph and addSeries_ . I have already read the ArcGIS help.

Can you guide me to some more examples and documentation?


texts:
I have been able to open WMS/WMTS layers normally until I found this one:
https://mapire.eu/en/map/thirdsurvey75000/wmts/
I don't understand why but it seems to be connected and when I add it to my map the title is shown at the side bar but the map won't load. I tried zooming in and out because some maps show up then. But with this one it just stays blank.


texts: FME: How can I use File Copy Writer to periodically copy new data within my project folders to a server?

texts:
All I want to do is copy a set of folders (and sub-folders) under a 'Projects' master folder on my laptop to the server. I want to automate this so it runs every hour or two and only copies new files to the server (so no overwrite).
Under the 'Projects' folder I have named project folders, i.e. 'Project A', 'Project B' etc.
The folder structure for each project is identical so they each have 'Incoming', 'Outgoing', 'Working' etc. and several sub-folders folders. Some are empty and some contain data. As project evolve, the standard set of folders gets populated with data and workflows.
The problem I have is that I can't seem to get file copy writer to copy this structure correctly. Data gets duplicated and the folder structure is lost. It also duplicates file gdbs at the root level of the destination directory and copyies ou the contents of gdb folders. Additionally it moves sub-folders to the root level so I end up with a mess...and duplicate data!
Should I run a workspace caller on each folder? Don't want a two step process so was hoping there was a way around using that approach.

[


texts: Formula to calculate optimal cell width and height of fishnet in ArcGIS?

texts:
I am trying to create a fishnet-layer in ArcGIS over my data and the program is asking me for the cell size. Now I know I can just type any data in there but I was told there is a mathematical formula to calculate the optimal cell size based on the area that I would like to cover.

I have tried to find it on Google but with no success so far, can anybody help me?


texts: Identifying gaps within polygons, and excluding islands in ArcMap

texts:
I have a digitised map of Ireland, consisting of 40 features where each feature is a different polygon representing an Irish constituency.

I have created a topology to help me identify gaps and overlaps in the map. The overlaps works fine, but unfortunately, Ireland consists of more than 800 islands (ironically). This means it comes up with a lot of gap errors (typically around 870), since these island polygons are not surrounded by another polygons.

The gaps I do want to find and correct are gaps within the mainland, between polygons, and donut holes. This means I have to go through 800+ errors individually to possibly find an actual error I want to correct.

I was wondering if there was another way of doing this, so I can ignore the island polygons, and find real gaps and errors within the mainland? Trudging through nearly 900 errors to find one actual error is tough work!


texts: How to update postgis layer from shapefile in QGIS?

texts:
My usual workflow with updating layer from shapefile was to:

Open DB Manager;

Load SHP files into QGIS;
Open DB Manager;
Import file (with drop option).

Now my problem is that I have created views on data in Postgresql and I can't just drop that table any more. Now my question would be:

If I can't drop tables any more how to effectively update PostGIS layer from SHP file ?

Since this is a point layer I currently make select by location / invert to get points not included in PostGIS layer.

Is there any effective solution to automate this proces or even better is there graphical tool to update PostGIS layer from SHP maybe based on attribute join / location checking?

Thanks !


texts:
I have 4 channel rasters in GeoTIFF.

I need to find the average value of these rasters and present them in the form of another raster, also in RGB, and I get a black and white file.

How can I get rid of this?
 gdal_calc.py -A /var/www/html/algorithm/cut103.tif -B /var/www/html/algorithm/cut1031.tif  --outfile=result.tif --calc="(A+B)/2"[![Вот что получается.][1]][1]

I also tried to use gdal_translate -of VRT input.tif input.vrt
My input.vrt:
<ColorInterp>Palette</ColorInterp>
<ColorTable>
    <Entry c1="170" c2="170" c3="255"/>
    <Entry c1="85" c2="170" c3="255"/>
    <Entry c1="66" c2="33" c3="18"/>
    <Entry c1="159" c2="81" c3="42"/>
    <Entry c1="205" c2="169" c3="21"/>
    <Entry c1="253" c2="254" c3="3"/>
    <Entry c1="230" c2="236" c3="6"/>
    <Entry c1="208" c2="223" c3="0"/>
    <Entry c1="185" c2="207" c3="2"/>
    <Entry c1="162" c2="192" c3="0"/>
    <Entry c1="138" c2="175" c3="0"/>
    <Entry c1="114" c2="160" c3="0"/>
    <Entry c1="91" c2="142" c3="3"/>
    <Entry c1="69" c2="129" c3="0"/>
    <Entry c1="45" c2="112" c3="0"/>
    <Entry c1="37" c2="96" c3="45"/>
    <Entry c1="21" c2="84" c3="45"/>
    <Entry c1="21" c2="68" c3="45"/>

Plese? someone, help


texts: Generating points along perimeter of circular feature in ArcMap?

texts:
How do I generate equally spaced points on the circumference of a circular feature in ArcMap?

I am currently working on a project where I need 36 equally spaced UTM coordinates that form a silo.


texts: Accessing shapefiles using Java in OpenShift

texts:
I have shapefiles in my resources directory and see them in the jar. After executing my service, the logger in OpenShift shows this: "Read lock: file:/opt/app-root/src/file:/opt/app-root/src/mars-de-service.jar!/BOOT-INF/classes!/CountryBoundaries.shp)". This tells me that it is able to find it and read it. In fact, in my local environment, the file is actually processed. However, in OpenShift, the next line in the logger shows this: "java.lang.RuntimeException: java.io.FileNotFoundException: /opt/app-root/src/file:/opt/app-root/src/mars-de-service.jar!/BOOT-INF/classes!/CountryBoundaries.shp (No such file or directory)". In my local, as I mentioned, I don't get such an error and my service works fine.
Below is my code:
URL res = CountryCodeLookup.class.getClassLoader().getResource("CountryBoundaries.shp");
File sourceFile = new File(res.getFile());
store = FileDataStoreFinder.getDataStore(sourceFile);
featureSource = store.getFeatureSource();
GeometryDescriptor geomDesc = featureSource.getSchema().getGeometryDescriptor();
attrName = geomDesc.getLocalName();

What am I missing?


texts: Copy Rows to SQL Server Table using ModelBuilder?

texts:
I want to copy a table in my GDB to a table in MSSQL server 2012.

I'm using the "Copy Rows" tool in ModelBuilder, but when I run it, I get the following error:

ERROR 000210: Cannot create output Database Connections\Connection to
<Server>.sde\<DB>.dbo.<Table>
DBMS table not found
DBMS table not found
Failed to execute

The Table exists, and I can even select it as the Output Table in the Copy Rows tool GUI.


texts: IDW Interpolation in QGIS gives weird values

texts:
I'm trying to do an IDW interpolation with QGIS tool.
But i'm really surprised by results. Interpolation values aren't similar as shape values.
I used a power of 0.2 to get a smooth raster.
Somebody have an idea ?

>>Here are data<<


texts: Matching two rasters with different projections and resolution

texts:
I have one elevation file (SRTM) and one satellite imagery file (sentinel-2), both are in GeoTIFF format.
The elevation file has a much larger area covered.
Here is the gdalinfo of the elavation file:

> Driver: GTiff/GeoTIFF Files: C:\DEV\3D\Assets\Heightmaps\nv\full.tif
> Size is 10801, 7201 Coordinate System is: GEOGCS["WGS 84",
>     DATUM["WGS_1984",
>         SPHEROID["WGS 84",6378137,298.257223563,
>             AUTHORITY["EPSG","7030"]],
>         AUTHORITY["EPSG","6326"]],
>     PRIMEM["Greenwich",0],
>     UNIT["degree",0.0174532925199433],
>     AUTHORITY["EPSG","4326"]] Origin = (-2.000138888888889,36.000138888888884)
> Pixel Size =(0.000277777777778,-0.000277777777778)
> Metadata:   AREA_OR_POINT=Area
> Image Structure Metadata:   INTERLEAVE=BAND Corner Coordinates:
> Upper Left  (  -2.0001389,  36.0001389) (  2d 0' 0.50"W, 36d 0' 0.50"N)
> Lower Left  (  -2.0001389,  33.9998611) (  2d 0' 0.50"W,33d59'59.50"N)
> Upper Right (   1.0001389,  36.0001389) (  1d 0'0.50"E, 36d 0' 0.50"N)
> Lower Right (   1.0001389,  33.9998611) (  1d 0' 0.50"E,33d59'59.50"N)
> Center(  -0.5000000,  35.0000000) ( 0d30' 0.00"W, 35d 0' 0.00"N)
> Band 1 Block=10801x1 Type=Int16,ColorInterp=Gray

and the satellite image:

>   Driver: GTiff/GeoTIFF Files:
> C:\DEV\3D\Assets\Imagery\new\sub_big_merge.tif Size is 15017, 8131
> Coordinate System is: PROJCS["WGS 84 / UTM zone 30N",
>     GEOGCS["WGS 84",
>         DATUM["WGS_1984",
>             SPHEROID["WGS 84",6378137,298.257223563,
>                 AUTHORITY["EPSG","7030"]],
>             AUTHORITY["EPSG","6326"]],
>         PRIMEM["Greenwich",0,
>             AUTHORITY["EPSG","8901"]],
>         UNIT["degree",0.0174532925199433,
>             AUTHORITY["EPSG","9122"]],
>         AUTHORITY["EPSG","4326"]],
>     PROJECTION["Transverse_Mercator"],
>     PARAMETER["latitude_of_origin",0],
>     PARAMETER["central_meridian",-3],
>     PARAMETER["scale_factor",0.9996],
>     PARAMETER["false_easting",500000],
>     PARAMETER["false_northing",0],
>     UNIT["metre",1,
>         AUTHORITY["EPSG","9001"]],
>     AXIS["Easting",EAST],
>     AXIS["Northing",NORTH],
>     AUTHORITY["EPSG","32630"]]
> Origin = (614090.000000000000000,3987110.000000000000000)
> Pixel Size =(10.000000000000000,-10.000000000000000)
> Metadata:
> AREA_OR_POINT=Area Image Structure Metadata:   INTERLEAVE=PIXEL Corner
> Coordinates:
> Upper Left  (  614090.000, 3987110.000) (1d44' 1.72"W,36d 1'18.55"N)
> Lower Left  (  614090.000, 3905800.000) (1d44'43.19"W, 35d17'19.93"N)
> Upper Right (  764260.000, 3987110.000) (0d 4' 6.98"W, 35d59'33.60"N)
> Lower Right (  764260.000, 3905800.000) (0d 5'42.84"W, 35d15'37.76"N)
> Center      (  689175.000, 3946455.000) (0d54'37.94"W, 35d38'37.76"N)
> Band 1 Block=15017x1 Type=Byte,
> ColorInterp=Red Band 2 Block=15017x1 Type=Byte, ColorInterp=Green Band
> 3 Block=15017x1 Type=Byte, ColorInterp=Blue

Now, how do I get a subset of the elevation file that exactly matches the same area of the satellite image with gdal_translate or any other tool (since as you can see the two files have differents projections, resolutions, units ... etc)?

In the meantime I have tried these two "solutions"

  Use GDAL command line to copy projections

  Use gdalsrsinfo to get the srs of the tiff that still has the
  projection:

  gdalsrsinfo -o wkt tiffwithsrs.tiff

  Then copy the output and use gdal_translate to apply it to a new tiff:

  gdal_translate -a_srs '...' tiffwithoutsrs.tif newfixedtif.tif

  just paste your projection after the -a_srs option
  wich did nothing .

and this which completely broke my file.


texts: Is it possible to reverse the trend/values of a raster?

texts: For example - if I have a raster DEM covering a stream with higher elevations in the west (up river) and lower elevations in the east (down river) is it possible to reverse this trend so that the east will have higher elevations? The idea is to essentially reverse the direction of the river.


texts: Managing sizes of input rasters for raster calculation in ArcGIS Desktop?

texts:
To calculate any formula (conditional statements/relational operation etc) in the raster calculator in ArcGIS using two different sources of raster images, I'd like to know input rasters must be the same size (row, col num) or not.

Generally, cell sizes of input rasters are the same but col and row numbers are slightly different.

Can the raster calculator operation manage their sizes automatically?


texts: Why does Calculate Geometry give incorrect number of decimals for Longitude?

texts:
intermediate ArcMap user here. I'm trying to calculate the geometry for a shapefile such that I can add lat/long columns to my attribute table. While I have no trouble calculating the latitude, the longitude keeps coming out a decimal place off.

I know that the accurate longitude for this shapefile should be -74.067, -74.068, etc. I'm not sure if the negative sign is screwing up the significant digits? I looked around StackExchnage and no one seems to have this problem.

For reference, when calculating the Latitude, for "Type", I used "float" with precision = 6 and scale = 5. This worked fine for Latitude but still gives me the incorrect result for Longitude.

Also for reference:
Projected Coordinate System:    PCS_CarMAGBOG
Projection: Transverse_Mercator
False_Easting:  92334.87900000
False_Northing: 109320.96500000
Central_Meridian:   -74.14659167
Scale_Factor:   1.00000000
Latitude_Of_Origin: 4.68048611
Linear Unit:    Meter

Has anyone faced this problem before, and if so, how did you resolve it? Is the problem in calculate geometry or in the projection itself?


texts: Appropriate data structure for a bunch of polygons

texts:
I'm not a GIS guy so this might be a very misguided question. At present I have compiled a python dictionary, on which I perform some simple processing and return some answers. My data is structured as a python dictionary right now:

entire_data = { "row1": {
                          "metadata1": "some info"
                          "metadata2": "some info"
                          "geometry": [(x1,y1),(x2,y2),...] #x,y are WGS84 coordinates
                        }
              }

The entire thing is a 100kb right now, but is expected to scale a lot more in the future. I plan to use this data (mostly) in the backend of a web application. Most of my toying around with the concept is done and now I need to somehow have it in a format for persistent storage. I looked into shapefiles but they appear like they might be overkill (also read some criticisms on the format). I can put them in a sql database, but the geometry field would have to be stored as varchar and that would mean a processing overhead each time an entry is to be fetched. What would be an appropriate format to put the aforementioned data into?


texts: Correlation analysis of two polygon shapefiles with nominal data

texts: I am working on a research- trying to model vegetation types based on factors affecting their distribution. In that case I have observed a strong relation between vegetation and soil types, from just overlaying the vegetation type polygon shapefile with the soil type polygon shapefile, now the problem is i want to quantify this relationship between the vegetation and soil types i.e.the vegetation-soil type correlation analysis so as to scientifically present this relationship. Please i want to know how to do that using ArcGIS 10.1, I have tried to search for solutions to this problem but with no success. I do not have the primary data which was used to create those polygon shapefiles so I am hoping there is way of just working with just the polygons.


texts: Post and Reconcile buttons in versionig toolbar grayed out in Arcmap

texts:
I have a versioned GDB (Postgresql) and I have versions DEFAULT-> A -> B.
Version A modifies two poly lines and version B modifies one of the two that is modified in version A.
In ArcMap, even though I have "List by source" selected and have selected the source (which also shows the selected version), I do not get "Post" nor "Reconcile" buttons enabled in the versioning toolbar when either version A or B is selected. Am I missing something or this is a bug?

For reference, I am using ArcGIS Desktop 10.3 with all patches installed as of now.


texts: Finding distance of a building from road junction

texts:
I have a peculiar problem which i cannot think of way to solve.

I need to find distance of a building from the nearest road junction. This would include distance to the nearest road and how far is that point of road from nearest junction.

I tried 'Near' but it gives either distance to road or distance to junction (displacement).

See image for example.


texts: Assessing impact of tall building on urban skyline?

texts:
I have a DSM of a small city center, and for planning purposes I wish to assess the impact of a very high building on the urban skyline. By impact I mean: given a tall building, what part of the urban landscape would not be visible from given viewpoints? While I have experience with viewshed analysis in ArcGIS (10.1), I do not manage to wrap my head around how to proceed on that specific issue. What I believe could be done is to select a number of (random) points within the urban area of interest, and calculate related viewsheds, in order to assess what part of the urban 'landscape' is actually blocked by the (future) high building. BUT, my issue is:

how can I simulate the presence of a high building (or, for instance, of added third or fourth floors on an existing building)?

Is it something that I can add to the DSM (how?)?

Can I use a shapefile, and, if I can, how the procedure should be set up?


texts: Simultaneously editing 2 features from 2 different layers in QGIS

texts:
Is there a way in QGIS to edit 2 different features from 2 different layers at once.

For example I selected a point in a point layer and a polygon in a polygon layer and I want to move them to another place at the same time.

I can enable editing on both layer but I still have to move my features one after the other which means I can't graphically move my 2 features by exactly the same x and y. That is a problem since I want them to stay at the same place relative to each other.

I found this post from 7 years ago saying it is not possible. I'm wondering if the situation has changed


texts: Storing features in PostGIS and then creating layers using queries

texts:
I have a whole set of line segments or tracks that represent real tracks or, in some case, boundaries between areas.

These segments form a network of lines on the map.

I have a number of GIS layers, e.g. walking tracks, vehicle tracks and a variety of zones (areas delimited by a set of segments) with a many to many relationship between segments and layers.

I figure the best way to deal with this is to put all the segments into a postgis table as individual features and then create a link table for each layer that lists then segments needed for that layer.

That way I have just one copy of each segment to maintain.

Is there a better way of doing this?
Can anyone point me to some tutorial material in this area.

More clarification:

Indeed, views (or sql queries - which amount to the same thing) were what I was assuming that I would use to build the layers.  I am pretty sure I know how to do what I want.

My questions was more along the lines of "Is the approach sound?".  As a GIS novice but with a extensive IT background  it seems to me that this must be a common situation i.e. having features that are used in multiple layers where you only want to maintain one copy of the feature.  I was wondering if the obvious approach (from a relational database point of view) was what people used in practice or if there were other approaches.


texts: qgis2to3 to migrate plugin from version 2 to version 3

texts:
I have a plugin in QGIS 2  and I'm trying to migrate it to version 3

I downloaded the 'qgis2to3' tool

--> I installed the pip (get-pip.py) :

python get-pip.py

--> I installed qgis2to3 :

python pip install --user qgis2to3

but when I try to use the command qgis2to3 /path/to/my/plugin it retruns this error :

can't find '__main__' module in 'qgis2to3'


texts: Merging lines connected by vertices (avoide connecting overlapping lines)

texts:
I have a layer representing hydraulic network and I would like to merge the lines when they are connected to each other at their start/end vertex.
In the picture below , I merged the lines manually by selecting the lines and using the tool 'Merge selected feature'. I would like to generalise the process to integrate it in a model builder.

I tried using 'Dissolve' then 'Multipart to singlepart' but this is not working. When lines have multiple branches, each of them becomes a new feature after splitting to singlepart (as in the picture in this thread)
I also tried using a 'Buffer' (with dissolve option activated) and then 'Joining by location' the lines within the buffer. This work for the lines in blue in the picture, however it doesn't work for the green and orange one because the dissolved buffer is joined where the green and orange lines are overlapping.
I tried to create a dissolved buffer that would except when lines overlap but no success.
I am working with QGIS 3.16.3


texts: ArcGIS JS-API click on Graphic and exceute code

texts:
Currently I can click on the GraphicsLayer and I can see an array of all my Graphics added there.

const layerMarkers = new GraphicsLayer({id: layerId});
this.map.addLayer(layerMarkers);
const graphic = new Graphic(point, imageSymbol, null, null);
layerMarkers.add(graphic);

But how can extract the data from that Graphic and do something with that data? I want to create a InfoTemplate only for that graphic.

layerMarkers.graphics[0].setInfoTemplate(infoTemplate);

I don't wan't to call .graphics[0] that way, I want to pick the array from clickin on the Graphic.

I was trying some this.map.on('click', (evt) => but what I need is to click on a specific Graphic.


texts: How to plot NDVI time series for two features on same plot?

texts:
I used the code mentioned below on Google Earth enfor plotting this. It shows error! Can anyone help with this. I want to plot NDVI time series for two features on same graph but two different lines.
var filteredCollection = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
  .filterBounds(table)
  .filterBounds(table2)
  .filterDate(start, finish)
  .sort('CLOUD_COVER', true);
var with_ndvi = filteredCollection.map(addNDVI);
var rgb_vis = {min: 0, max: 0.3, bands: ['B4', 'B3', 'B2']};
Map.addLayer(filteredCollection.median(), rgb_vis, 'RGB');
var greenest = with_ndvi.qualityMosaic('nd');
Map.addLayer(greenest, rgb_vis, 'RGB (greenest pixel)');

print(Chart.image.series(with_ndvi.select('nd'), table&table2));


texts: Clustering not working on features in OpenLayers

texts:
I've built a mapping solution using OpenLayers that works well except I can't get the points to cluster.  I've looked at the examples on the site, and I've even implemented those examples in my environment to make sure they work.
When I try to do it with my map, however, all of my points show up, but they're not clustered.
Here is my code where I create the Vector Layer:

var layerVector = new OpenLayers.Layer.Vector("Regions",{
                 strategies: [

                        new OpenLayers.Strategy.Cluster()
                    ]
            });

I add all of the points to an array of features:

var feature_point = new OpenLayers.Feature.Vector(
                            this, {
                                "rName" : rObj.tName
                            }, {
                                'fillColor' : rObj.rColor,
                                'strokeWidth' : borderWidth,
                                'strokeColor' : rObj.rColor,
                                'fillOpacity' : regOp,
                                'pointRadius' : regPercent
                            });
                            fArray.push(feature_point);

Add the features:

layerVector.addFeatures(fArray);

And add the layer to the map:

window[mapName].addLayer(layerVector);

Seems pretty simple, but it seems like I'm missing something.
In the example on the OpenLayers site, they add a styleMap to the layer, but I tried this and it doesn't do anything for me.  I'm honestly not sure what it even means:

 var style = new OpenLayers.Style({
                    pointRadius: "${radius}",
                    fillColor: "#ffcc66",
                    fillOpacity: 0.8,
                    strokeColor: "#cc6633",
                    strokeWidth: "${width}",
                    strokeOpacity: 0.8
                }, {
                    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;
                        }
                    }
                });

                strategy = new OpenLayers.Strategy.Cluster();

                clusters = new OpenLayers.Layer.Vector("Clusters", {
                    strategies: [strategy],
                    styleMap: new OpenLayers.StyleMap({
                        "default": style,
                        "select": {
                            fillColor: "#8aeeef",
                            strokeColor: "#32a8a9"
                        }
                    })
                });

Any help would be immensely appreciated.


texts: Passing user created x,y coordinates as input feature to identity task using ModelBuilder?

texts:
Is there a way to pass user created x,y coordinates as input feature to the identity task?

I don't want to choose a table with x,y fields nor a feature class.


texts: Disconnect QGIS from database without closing the whole program

texts: I'm making PostGIS queries to a GIS database from python and I'm using QGIS to visually check the results of my queries, but every time I want to run a new query I have to close QGIS because it keeps its database connection open. Do you have any idea how to force QGIS to disconnect from the database without closing it every time?


texts:
Trying to load the service http://gis.ktimanet.gr/wms/ktbasemap/default.aspx
and receives the message
Could not understand the response. The wms provider said:
Could not get WMS capabilities: error occurred while parsing reference at line 5 column 179
This is probably due to an incorrect WMS Server URL.
Neither in http://gis.ktimanet.gr/wms/ktbasemap/default.aspx?SERVICE=WMS&REQUEST=GetCapabilities
I have result.I never had a similar problem again with this service into something so simple


texts: Automatically converting GeoPackages to ZIP archives with Shapefiles and QML/SLD with native QGIS

texts: Currently we use a third party tool (mygeodata website) to convert our geopackage back to a zip folder with all the contents of the geopackage and the styles intact (even the standard Style Setting). But we want this to do purely in QGIS and not just a simple export but all of the geopackage being exported at once. Is this possible?


texts: Shifting z values of all features in a shapefile?

texts:
I have a few line and polygon shapefiles that have the wrong z vertex values. I need to apply the same z shift for all features in all the shapefiles.

I looked around and only managed to do it by loading the data in a PostGIS database and running the following SQL script on every layer:

UPDATE public.poly
SET geom = ST_Translate(geom, 0, 0, 10)

Is there an easy way to shift z values of all features in a shapefile in QGIS, even if it involves using PyQGIS?


texts: How to overlay map layers (shp and csv) in R?

texts:
This is a crosspost from stackoverflow. I know very little about GIS software and am doing all my mapping in R. Apologies in advance if this is too basic question. Let's say I have two shapefiles from different sources but with different attributes. Say, one is for adminstrative boundaries of Texas (boundaries.shp) and the other is for Texas rivers (rivers.shp). I also have a third file towns.csv showing the locations of towns within the state. After reading in the files, I can overplot town locations on the adimistrative boundaries in the maptools package:

plot(boundaries); points(towns$lon, towns$lat)

But how can I overlay all three? Surely there is an easy way of doing this?


texts: How to overlay Radiance (Top-of Atmospheric radiance) over Digital Numbers Image?

texts: I am trying to process INSAT-3D images, which contain six different bands. I have to extract each and every band form the .h5 data file, so when I extract first band (MIR) band form the data set, I found that My MIR band image size is (3062, 3207) which digital number ranges form 0 to 1023, while its radiance size is (1, 1024) which ranges form 0 to 0.3. Can anyone please tell me how could I overlay these radiance values onto my image.


texts: Add Custom Data to OpenLayers Plugin for QGIS

texts:
Is there a way to add custom data to the OpenLayers plugin for QGIS?  I have tiled data using Map Tiler, which has produced my zoom levels, and an OpenLayers.html.  I am hosting it internally on our servers, and I would like to access it from the OpenLayers plugin.

Thoughts on where to start?

Thanks...


texts: Arcgis Server 10.1: Workflow/Setup for Secure Content?

texts:
We have a server which has ports 80 and 6080 open to the outside world--and probably port 443 as well. But port 6643--the default secure port for ArcGIS Server is not open to the outside. The server also has a functioning Web Adapter for port 80 based calls to Services. This is a Windows 2008 Server.

My supervisor wanted me to install an SSL certificate on urgent basis; so I have enabled both HTTPS/HTTP on the server and also created and sent a .csr file to our IT staff so they get as a CA issued SSL certficate. But for now at least we can, inside the network, get to  even with see the security warning and browser exception option.

But, according to this: https://geonet.esri.com/thread/61639 we should not install the CA issued SSL on our server? If not then how to protect content? Let's say if there are applications accessing Services and logins to port 6080 or 80 then they remain not secure. So we will need to have a genuine SSL certificate in the picture. But how?

This is my first time dealing with ArcGIS Server's SSL layer and so my apologies if I am way off. Just trying to understand.


texts: Selecting features enclosed by another layer using QGIS?

texts:
I have 2 layers selected layers as shown in the figure below, hotosm_glp_rpads_lines (Layer 1) and LaSoufriere_Imapact_Zones(Layer 2). I want to select all the polylines in layer 1 which is bounded by layer 2 (all the polylines which are inside the brown region).

How do I do this in QGIS?

NOTE: It would be great if there is some code that can be written in the select feature by expression under the attribute table since I have to also select another set of polylines. In more detail, I need to,
I need to select polylines bounded by the brown region & also select polylines having the attribute "highway" = "primary" in the attribute table


texts:
I want to only put labels on the vertices of the polylines.  I can use the data defined features to place a label at specific points, but not at multiple specific points.

Is it possible to label only the polyline vertices in QGIS?


texts:
I am blocked on trying to work with Layer Properties Rule-based Expressions.

For a simple example in expression string builder:
I have a column field named HOUSE.
A list of values within that field show SMALL, MEDIUM, LARGE.
What rule expression would I use to allow different colors for each of the values?

I tried this:

CASE WHEN "HOUSE" = 'SMALL' THEN color_  rgb(255,0,0) END

but got this message:

  Parser Error: syntax error, unexpected COLUMN_REF, expecting WHEN or
  ELSE or END Eval Error: No root node! Parsing failed?

Can someone please reply with a correct expression I can use?


texts: Is there a way to get offline tile data for Cesium?

texts:
I would like to work with Cesium (www.cesiumjs.org) without any internet connection. The documentation provide a guide for working offline and prvides a way to import tiles (see  https://github.com/AnalyticalGraphicsInc/cesium/wiki/Offline-Guide). However this way does not import data to be used by the TerrainProvider.

Is there a way to provide the data locally or via local cache?
I know of the small terrain from http://cesiumjs.org/data-and-assets/terrain/small-terrain.html , However I need a finer terrain detail.


texts: Problem with photo displaying in QGIS 3.10

texts:
I would like to know how to display a photo in a image viewer with a featured action.
I create a layer with 3 fields: id, leyenda and filepath. Filepath is the field where I put my photo..

Then I created an action to open the picture who will be attached

But when I use the action (the cursor change) nothing happened.
Do you know why and how to solve the problem ?


texts: Displaying geometry returned from ST_Buffer in PostGIS with OpenLayers?

texts:
I'm using PostGIS, GeoServer and OpenLayers, and the thing is that I want to display the geometry result from a query like this:

SELECT ST_Buffer(ST_GeomFromText('POINT(100 100)'), 500);

To perform queries I'm using Java, the line I use to perform a query is the following:

resulset = statement.executeQuery(query);

where resulset is a ResultSet type object, when I convert it to String I get the following:

org.postgresql.jdbc4.Jdbc4ResultSet@163b102

My question is how do I interpret that to display the POLYGON I suppose is in there in a layer? or am I doing something wrong? Hope someone can help me with this question.


texts: How to store a vector layer as a geojson file in openlayers 3?

texts:
I have some vector layers that I import from geojson files and run some search on them, I want to display only the results of the search as a vector layer.

I am using a jquery form to search, and I get the result but it doesn't show on the map because the page is after that refreshed and the result is lost. I thought of exporting the result in a geojson file this way I can preserve it.

Here is the code for the search it works fine except that the result layer isn't added to the map.

jQuery('#default-search').submit(function(event) {

    var Results=[];
    var src=layer.getSource();
    src.forEachFeature(function(feature) {
        if(feature.get('Name')==='some_Name'){
            Results.push(feature);
        }
    });
    tmp_source=new ol.source.Vector({
                features: Results
                });
    tmp=new ol.layer.Vector({
            title:'temporary layer',
            source: tmp_source
            });

    map.addLayer(tmp);
    tmp.setVisible(true);

});

Is this possible or is there another way easier ?


texts: Finding the Center Line from a set of 3D Points

texts: I have a set of 3D points. They follow a curved pattern with a rather constant diameter as shown below. What would be the algorithm to trace the approximate center line of these points?


texts: How to decrease GeoTIFF file size by eliminating or negating the blackfill?

texts:
I have a 5 band, 16-bit orthorectified GeoTIFF from a satellite image that contains a significant amount of blackfill and has a file size of over 21GB. I would like to decrease the file size by eliminating or negating the blackfill. I have attempted to use the SPARSE_OK option in GDAL as such:

gdal_translate -of GTiff -co "BIGTIFF=YES" -co "SPARSE_OK=TRUE" -co "TILED=YES" input.tif output.tif

This has had no effect on the file size. Is my understanding of how SPARSE_OK works incorrect or is there a potential problem with the SPARSE_OK option? Is there a better way of reducing the file size? I have tried "COMPRESS=DEFLATE", which works, but then the product is veryyyy slow is loading into the viewing software.


texts:
The numeric scalebar in the map composer shows too detailed information on the output.
Moving the Map canvas changes the scale depending on the moved distance on the map.

We use EPSG 4258 (ETRS89)

How do I round the shown output by -2 digits? It should round to hundreds


texts: PostGIS Out of db VRT containing JPEG2000 GDAL 2.1 OpenJPEG

texts:
Build GDAL 2.1 with OpenJPEG
Build PostGIS with above GDAL build
Created VRT file with SENTINEL2 Driver / JP2 files.

gdal_translate SENTINEL2_L2A:S2A_USER_MTD_SAFL2A_PDMC_20160520T042556_R030_V20160520T000224_20160520T000224.xml:10m:EPSG_32756 10m.vrt -of vrt

raster2posgresql the VRT into PostGIS with -R flag (out of raster db using full path to file)

raster2pgsql -R /home/ben/S2A_USER_PRD_MSIL2A_PDMC_20160520T042556_R030_V20160520T000224_20160520T000224.SAFE/10m.vrt -b 2 -C -t 128x128 -P test.T56JKN_20160520_A_10m_b02 | psql -d giddb

When I query:

select
     ST_Value(rast, ST_Transform(ST_PointFromText('POINT(149.9878 -29.3507)', 4326),32756))
from
     test.T56JKN_20160520_A_10m_b02
where
     ST_Intersects(rast,ST_Transform(ST_PointFromText('POINT(149.9878 -29.3507)', 4326),32756))

I get error:

  ERROR: rt_raster_from_gdal_dataset: Could not get data from GDAL
  raster SQL state: XX000 Context: PL/pgSQL function
  st_value(raster,integer,geometry,boolean) line 18 at RETURN

Other notes:
Out of db works perfectly fine with a VRT containg a tif file and also with a tif file on it's own.

Any advice or comments welcome. I am not sure if this is even possible?

Edit: Would using the Kakadu JPEG2000 driver help?

UPDATE:

OK so I've built everything from scratch again and can confirm:

JP2OpenJPEG is listed in PostGIS query: "SELECT short_name FROM ST_GDALDrivers();"

JP2OpenJPEG is available in GDAL by running "gdal_translate --formats"

I can gdal_translate a SENTINEL2 driver created VRT containing .jp2 to GTiff

PostGIS version information:
POSTGIS="2.2.2 r14797" GEOS="3.5.0-CAPI-1.9.0 r4084" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 2.1.0, released 2016/04/25" LIBXML="2.9.1" LIBJSON="0.11.99" TOPOLOGY RASTER

When trying to query raster I now get error:
ERROR:  rt_raster_from_gdal_dataset: Could not get data from GDAL raster
CONTEXT:  PL/pgSQL function st_value(raster,integer,geometry,boolean) line 18 at RETURN

Other notes:
When I run raster2pgdsql with -s [dataset epsg] I get error:
ERROR:  rt_raster_from_gdal_dataset: Could not get data from GDAL raster
CONTEXT:  SQL function "st_transform" statement 1

When I run without -s flag it completes OK with and without tiling.


texts:
We work with lots of DXF files, and it takes much time to save imported DXF file as SHP or convert them to polygons (even using batch converter) to enable "Edit" button.

Is there any way for quicker editing lots of DXF layers right after import?

And quicker way to make qgis recognize them as polygons?


texts:
I'm using ArcGIS 10.0 - ArcInfo license. I'm a fairly new user.

I'd like to align rectangular symbols perpendicular to to the lines that reside nearby.

I have two layers. A layer containing lines and a layer containing points. Both are shapefiles.

In this case I'm using NHDFlowlines (streams and rivers) and the rectangular symbols will represent dams.

I'm attempting to use the 'Align Marker to Stroke or Fill' tool.

I've taken each shapefile and exported them out to Personal Geodatabases.
I've applied which symbols I want to use for the layers.
I then use "Convert Symbology to Representations".

Here's where I think I have my first problem: When the new representation layer gets added, the check box is grayed out and there is a little red exclamation point next to it.

My next step was to use "Repair Data Source" where I point to the file within the geodatabase.
I then save my two Representation layers as .lyr files.
After that, I try to run the 'Align...' tool.

It processes successfully, but nothing actually aligns...

I've gotten to this point through trial and error, but I'm not sure what doing...


texts:
I want to use controls like DrawFeature and SelectFeature, but while these controls are active, I want to have a way to drag the map, using DragPan control.

For example, I am selecting features. I select then first, the I would have to pan the map in order to select the 2nd feature (not unselecting the 1st one).

Is it possible to make DragPan work as some kind of background control, that would always be active, even when using other controls?


texts: How to split and extract multispectral imagery in QGIS (2.18)

texts:
I am looking for a way to split multiband images and extract all the bands at once- instead of subsetting one band at a time.

Is there any method in QGIS which could allow for this? I know there are split RGB band tools but I have 8 band files.


texts:
I may be facing a unique issue. My ArcGIS server is running within my company network.The server has picked up the proxy settings of the box it is running in.

But to serve the printing request, the server bypasses this proxy settings and calls url(services.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer) directly. As a result this request is getting blocked at the company firewalls and the request is failing. When the firewall was disabled the printing is a success.

I have followed this link
http://resources.arcgis.com/en/help/main/10.1/index.html#//015400000309000000

How do I route all the server requests to the Internet through my proxy server?

I am running ArcGIS 10.1 Server on Windows.


texts: How to convert FEMA flood map (FEMA NFHL v3.2 kmz) to GeoJSON

texts:
TLDR: Using Python, how can I convert a kmz file to GeoJSON? It appears that the kmz file makes web service calls.

Details:

I downloaded the FEMA_NFHL_v3.2.kmz file from FEMA here:
https://hazards.fema.gov/femaportal/wps/portal/NFHLWMSkmzdownload

My goal is to take that file and convert it to a format that I can feed to MapBox or otherwise use to populate a PostGIS database. My preferred format is GeoJSON.

I found the following code on the StackExchange network, but it doesn't seem to work.

import gdal

srcDS = gdal.OpenEx('./FEMA_NFHL_v3.2.kmz')
ds = gdal.VectorTranslate('output.json', srcDS, format='GeoJSON')

The code above produces this error:

$ python scratch.py
Warning 1: The output driver does not natively support DateTime type for field timestamp. Misconversion can happen. -mapFieldType can be used to control field type conversion.
Warning 1: The output driver does not natively support DateTime type for field begin. Misconversion can happen. -mapFieldType can be used to control field type conversion.
Warning 1: The output driver does not natively support DateTime type for field end. Misconversion can happen. -mapFieldType can be used to control field type conversion.
ERROR 1: Layer 'NFHL' does not already exist in the output dataset, and cannot be created by the output driver.
ERROR 1: Terminating translation prematurely after failed
translation of layer NFHL (use -skipfailures to skip errors)

I unzipped the kmz file, and have the following:

$ tree FEMA_NFHL_v3.2
FEMA_NFHL_v3.2
├── doc.kml
└── files
    ├── FEMA_Logo.png
    ├── GenStruct_item.gif
    ├── LOMA_item.gif
    ├── Levee_item.gif
    ├── LiMWA_item.gif
    ├── NFHL_Legend.png
    ├── NFHL_Quick_Start_Instructions.png
    └── RivMrk_item.gif

1 directory, 9 files

It looks like the only file worth considering that contains any data is doc.kml. Inspecting doc.kml, I see there's some code referencing a FEMA webservice. Further, opening FEMA NFHL v3.2 kmz in Google Earth, I see that it is, indeed, making web service calls.

Is there a way to programatically, using Python, have the .kmz file make the webservice calls and then generate a GeoJSON file?


texts: Gathering details of Annotation Subclass from ArcGIS Pro map.ListLayers() objects

texts:
I am attempting to remove all Annotation Subclass layers from a list that is generated from a Map inside of ArcGIS Pro.
From what I can tell the subclasses must be a "Layer" inside of Pro.
I cannot find a method to conclusively distinguish which layers are annotation subclasses.
l = [list of Layer Objects]
for x in l:
    if x:
        try:
            desc = arcpy.Describe(x)
            print (desc.featureClass.featuretype)
        except:
            print (f"{x} cannot be described")

This will yield "Annotation" for both the Annotation layer, and its subclasses.
desc.name

This will yield the Name of the Annotation layer for all subclasses
for x in l:
    print (x)

generates an error NameError: The attribute 'longName' is not supported on this instance of Layer as soon as it hits the Annotation subclass


texts: Python in ArcGIS: How to add a space between upper and lowercase letters in TOC

texts:
I'd like to use python to update the TOC by adding a space (white space) for layer names that are like this: "BoundaryBuffer"  or "WaterEdge" and replace them with "Boundary Buffer" or "Water Edge".

I'm new to python and could use some assistance so far I've gotten to this but I can't get it to work correctly:

    import arcpy
>>> import re
>>> mxd = arcpy.mapping.MapDocument("CURRENT")
>>> df = arcpy.mapping.ListDataFrames(mxd)[0]
>>> layers = arcpy.mapping.ListLayers (mxd)
>>> for lyr in layers:
...     re.findall('[A-Z][a-z]*', lyr.name)
...     arcpy.RefreshTOC()

I've looked at this post as well but I don't know how to integrate it with lyr.name in ArcGIS.

EDIT: I solved it myself! If anyone is interested or wants to comment:

import arcpy
... import re
... mxd = arcpy.mapping.MapDocument("CURRENT")
... df = arcpy.mapping.ListDataFrames(mxd)[0]
... layers = arcpy.mapping.ListLayers (mxd)
... for lyr in layers:
...     tempName = re.sub(r'([a-z])([A-Z])', r'\1 \2', lyr.name)
...     lyr.name = tempName
...     arcpy.RefreshTOC()


texts: QGIS Spatial join many points to single vector

texts:
I've got two shape files, one contains many point coordinates and the other vector shapes. There are many point coordinates contained within each vector shape.

I would like to perform a spatial join and create a new shapefile of the coordinates but with the attributes of the overlying vector shape.

QGIS spatial join does not seem to work for this. Any ideas? Also setting up a spatial database would be the ideal solution but is not possible in the time frame required.

Thanks
Richard