texts: Selecting or highlighting upstream elements

texts:
I would like to use the information from the Swedish SMHI related with their hydrological model elements (catchment areas) that have a Containment hierarchy structure, each element has the information about its dowstream element, to select or highlight a whole catchment area like in the picture below:

The first step is to select them, but once that is done, of course there are other question.

Each element, among other data, has AROID (the ID of the element) and OMRID_NED (the ID of the element downstream).

I would like to do this in QGIS in a simple way, like a click in one element.

The data I want to use is here:
Delavrinningsområden SVAR_2016_3 (71,0 MB, zip)
An example of how they make it to work on their webpage is here SMHI vattenwebb Modelldata per område


texts: QGIS: Every input to QuickWKT results in a point at null island

texts:
If I open the QuickWKT plugin in QGIS, and choose "POINT (WKT)" from the combo box, and enter the text,

POINT (30 10)

I get a point at "null island" (0°N, 0°E). The entered WKT appears well-formed to me.

What's going wrong?

Note: this answer includes a screenshot of QuickWKT. My QuickWKT does not look like that; it looks like this:

It isn't clear to me as a user why I need to specify POLYGON / POINT / etc. when the WKT itself encodes that. But the examples it has have those keywords, so I presume I'm supposed to enter them as well.


texts:
I want to reproject data via sld:

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0"
 xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
 xmlns="http://www.opengis.net/sld"
 xmlns:ogc="http://www.opengis.net/ogc"
 xmlns:xlink="http://www.w3.org/1999/xlink"
 xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <NamedLayer>
    <Name>RAIN_SLD</Name>
    <UserStyle>
      <Title>RAIN</Title>
      <FeatureTypeStyle>
        <Transformation>
          <ogc:Function name="gs:BufferFeatureCollection">
            <ogc:Function name="parameter">
              <ogc:Literal>features</ogc:Literal>
              <ogc:Function name="gs:Reproject">
                <ogc:Function name="parameter">
                  <ogc:Literal>features</ogc:Literal>
                </ogc:Function> <!-- end parameter -->
                <ogc:Function name="parameter">
                  <ogc:Literal>forcedCRS</ogc:Literal>
                  <ogc:Literal>EPSG:4326</ogc:Literal>
                </ogc:Function> <!-- end parameter -->
                <ogc:Function name="parameter">
                  <ogc:Literal>targetCRS</ogc:Literal>
                  <ogc:Literal>EPSG:900913</ogc:Literal> <!-- end parameter -->
                </ogc:Function> <!-- end targetCRS -->
              </ogc:Function> <!-- end reporject -->
            </ogc:Function> <!-- end parameter -->
            <ogc:Function name="parameter">
              <ogc:Literal>distance</ogc:Literal>
              <ogc:Function name="env">
                <ogc:Literal>BUFFER_DIS</ogc:Literal>
                <ogc:Literal>10000</ogc:Literal>
              </ogc:Function>
            </ogc:Function> <!-- end parameter -->
          </ogc:Function>
        </Transformation>

        <Rule>
          <Name>BUFFER</Name>
          <Title>BUFFER</Title>
          <ogc:Filter>
            <ogc:PropertyIsEqualTo>
              <ogc:PropertyName>SN</ogc:PropertyName>
              <ogc:Literal>A0001</ogc:Literal>
            </ogc:PropertyIsEqualTo>
          </ogc:Filter>

          <PolygonSymbolizer>
            <Stroke>
              <CssParameter name="stroke">#d3d3d3</CssParameter>
              <CssParameter name="stroke-width">2</CssParameter>
            </Stroke>
          </PolygonSymbolizer>

          <PointSymbolizer>
            <Graphic>
              <Mark>
                <WellKnownName>circle</WellKnownName>
                <Fill>
                  <CssParameter name="fill">#ff0000</CssParameter>
                </Fill>
              </Mark>
              <Size>6</Size>
            </Graphic>
          </PointSymbolizer>
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>

When i requested the layer with the SLD above, i got an error:

  Could not convert the value EPSG:4326 into the expected type interface
  org.opengis.referencing.crs.CoordinateReferenceSystem for parameter
  forcedCRS.

I found an example:
Geoserver: Georectification via SLD fails

It is useful. But in my situation, 2 CRSs(forcedCRS and targetCRS) must be set.

How could I do with it?


texts: Displaying .PNG file folders generated from MBUtil in HTML Doc

texts:
I recently followed this tutorial http://blog.davidelner.com/create-map-with-tilemill-and-leaflet/ that illustrated how to generate PNG tiles from Mbtile data sources. I was successful in creating the PNG tile folder dictionary on my local machine with the MBUtil utility. I wanted to know how would I serve the tiles in my .HTML document for offline use. This is what I currently have in my script but have been unsuccesful in displaying the tiles.

  var newLayer = L.tileLayer('file:///C:/Users/067626/Documents/GitHub/mbutil/libya/{z}/{x}/{y}.png', { maxZoom: 8 });
  map.addLayer(newLayer);


texts: Retrieve the attributes from a feature with IdentifyTask/IdentifyRequest (Javascript API)

texts:
So I'm building a GIS web application.  I am able to get the result from IdentifyTask with the "oncomplete" function, the problem I am running into is how I can get the information back out of the feature/attribute.  Here is the code that I am currently working with:

function onMapClick(evt) {
    var identifyParams = new esri.tasks.IdentifyParameters();
    identifyParams.geometry = evt.mapPoint;
    identifyParams.tolerance = 2;
    identifyParams.mapExtent = _Map.extent;
    var identify = new esri.tasks.IdentifyTask(url);
    identify.execute(identifyParams);
    identify.on("complete", identifyCallBack);
    identify.on("error", showError);

}

function identifyCallBack(lotResults) {

    var out = [];
    var result = lotResults.results[0];
    var feature = result.feature;
    var attribute = feature.attributes;

    dojo.forEach(attribute, function(attr){
        out.push(attr);
    });
    alert(out.join(", "));

}

function showError(results) {

    alert(results.error.message);
}

Right now my identifyCallback is just coded trying to get the data I want.  Would it be easier to do a dojo.postback and pass the result back to my ASP.net form and parse up the data in Visual Studio?  If so how would I code the dojo.postback to return the value to a RadAjaxManager in Visual Studio?


texts: Legend symbol is not visable in Layout Composer when styling points with two different attributes

texts:
In QGIS 3.10.2 I have a map that I want to show points by colour using one category and size by another. I have tried to follow the instructions here: https://docs.qgis.org/3.10/en/docs/user_manual/working_with_vector/vector_properties.html#proportional-symbol-and-multivariate-analysis. This I have done but a consequence is that in the Layout composer the legend symbol is either missing or very tiny.
This is my map with just a colour classification:

I then add the instructions to re-size by a second attribute (I converted a 4 category variable into the numbers 1, 2, 3 and 4 and called it SizePoint)

Now the legend in the Layout composer does not show a dot at all

Now here this isn't a problem because all the dots are of the same INSTANCE. But there are other maps where there are more than one instance:

You will see that there is no indication in the legend what INSTANCE 1 or INSTANCE 2 are.
I did try and look at the very specific help text here:

But I do not know what the classification frame is or see a Change button.
Any suggestions?


texts: Copy protection via license file for ModelBuilder models?

texts:
I've build a large model in the ArcGIS ModelBuilder and I want to share this toolbox with some other people BUT I don´t want to lose control over it. The first step was to set a password to my script, but this only protects the script not the functionality. If someone have my toolbox, he could use it when and where he wants and share it with others like he wants. That's not the way I want it.

Is there a way to build something like my own license file? My idea was to use python (or any other option) at the beginning of the model to check if there is an license file on my Gdrive account, something like an *.txt or *.html. The model only should work, if this file was found (I´ve done that on a local machine - and it works - but not with using the cloud).


texts: Joining attributes from points to polygons

texts:
I have a polygon layer without data, and a point layer with attributes.

Each point is contained in one polygon, and I would like to have the data from the point layer joint in the polygon layer, so every polygon has the data from the point contained in it.

Is this possible with QGIS?


texts:
I'm using MapBasic and have used a query to get the rowids (populated in a column) of minimum values

Select SUBSECT_ID, Row, Min(MaxLV)from LVQuery group by SUBSECT_ID into MinLV

Now, I want to delete these entries from the main Table (LVQuery). However I cannot as MapInfo won't let me delete from a view. How do I circumvent this.


texts: Hiding / encrypting password in QGIS processing script

texts:
Is it possible to hide or encrypt password in QGIS processing script?
For example, hide password for PostgreSQL connection.
QGIS "Master password" is possible, but I need a solution for code of processing script.


texts:
I'm trying to mosaic two raster images using R. I used this script:

setwd("D:\\Modis EVI\\")
raster_dir <-    "D:\\Modis EVI\\"
r1 <- brick(paste(raster_dir,"evi8_stack.tif", sep=""))
r2 <- brick(paste(raster_dir,"evi9_stack.tif", sep=""))
m1 <- mosaic(r1, r2, fun=mean)

The script worked fine but the result mosaic has this black line, see picture

Any advice please on why I have this line, and how could it be removed?


texts:
I have managed to load QgsVectorLayer from PostgreSQL and add it on the map. Is possible to trap save event after editing, do my own saving, call layer refresh and in this way prevent from default save by data source?
The idea is to have reading through the standard provider and save changes through the application server applying certain logic and checks.


texts: Publishing to ArcGIS Server using ArcPy and AnalyzeForSD()?

texts:
I'm trying to publish some data layers to ArcGIS Server, using ArcPy.  This is my first time trying this.

I can create my feature classes, and save them in a file geodatabase.  I can then create layers corresponding to them.  I then save the layers to individual layer files, and thence to an mxd.

I then create an SD draft document (using CreateMapSDDraft).  No problems so far.

The next step seems to be calling AnalyzeForSD, passing in the SD draft.

Two things at this point.  One is that that is supposed to, as near as I can tell, return a Dictionary with 'errors', 'warnings', and 'messages'.  However, I was having some issues and was getting an IOError that had a Dictionary just like that attached to it.

Was I doing something wrong there (other than the problems that were causing the IOError)?

Second, once I fixed those problems, I started getting another IOError, but the explanation for this one was:

  Unidentified error


texts:
I have a point file with elevation data. Some of the points are elevations relating to soil horizons.

I've created a TIN surface for each soil horizon then used the "extrude between" tool to create a 3D polygon of the area between each surface, basically creating a set of blocks stacked on each other representing the different soil horizons.

I now want to select points that fall within each soil profile. However, my multipatch is not closed so using the tool "inside 3D" does not give me all the points that fall inside each soil profile. I have about 5 different profiles, one of which is closed and the others are not. I've tried the "enclose 3D" tool but it adds volume in the z direction, misrepresenting the soil profile.

Why is the extrude between tool not enclosing all multi-patches and how can I fix this?


texts: Publishing service using service definition file (sd) through the ArcGIS REST API?

texts:
The algorithm I have so far is basically this:

Create a sddraft from a mxd file (using arcpy.CreateMapSDDraft().
Stage sddraft into a sd file (using arcpy.StageService_server().
Upload sd file to ArcGIS server over REST.

There's a little dance that has to happen due to the fact that we use web-adaptor and have to redirect all requests to a single-sign-on server first for authentication, before forwarding the intended request to the intended destination.

Call createService on ArcGIS server's REST API, pointing to the uploaded sd file.

This happens in the same http session, so no need to go through the authentication dance again.

The fourth step manages to create the service, but startup fails due to Invalid file. The only documentation for the createService endpoint that I could find is this, and I can only surmise that the endpoint only accepts the old and deprecated msd file, but not the newer sd file. If that's the case, what endpoint am I supposed to hit for creating a service from a sd file? Note that arcpy.UploadServiceDefinition() isn't an option for us due to the fact that authentication is handled elsewhere from the ArcGIS server.


texts:
I'm planning to create an online map with openlayers (2 or 3) or leaflet. I have to show some points/markers (~30.000), polygons (~500) and lines in an area of ~10x10 kilometers in several layers. This data comes in local x and y coordinates. The background will be made of raster tiles, perhaps OSM.

Now I'm wondering, whats the best strategy to do this.

Is it possible to make the local coordinate system the default for the whole map in one of those map libraries?
If not, should the data be converted into the coordinate system of the map (WGS 84) in the database or a cache or doesn't that matter?


texts: Calculating area of residential area wihin multiple ring buffer in ArcMap?

texts:
I need to calculate the size of residential area within 0-250m, 250-500m and 500-750m of a bus stop.

I have made a multiple ring buffer around the bus stops. And I have a shapefile of the residential area.

But I can't work out how you can combine the shapefile of residential area with the multiple ring buffer and then calculate the size of residential area of each buffer ring.

Does someone know how to do this?


texts: Spatial Index on 100s of Millions of Records Taking Days in PostGIS?

texts:
I have a table with ~200million spatial records in it, and it is taking days to build a spatial index on the one geometry column: should this be taking so long?

The geom column is of the type geometry(Linestring) and is representing features on chromosomes.  I used copy to load these data in roughly 10 minutes, but when I try to build a spatial index on the geometry column, it appears to spin indefinitely (specifics on my setup below the questions). Building a standard index on a table representing the Linestring not as geometry but as three columns: y_coord, x_start, x_end only takes minutes... I feel like I'm doing something wrong in building the spatial index.

My questions are:

Should a spatial index on 200 million records take days to generate?
What is the best configuration to build spatial indexes on datasets of 100s of millions to billions of data points?

My table looks like:

CREATE table JEM2(
idx serial,
var geometry(Linestring),
record text);

I load the 200million rows into the table with copy. Records loaded into the table look correct:

SELECT * from jem2 where idx=1;
 idx |                                        var
          |                   record
-----+--------------------------------------------------------------------------
----------+---------------------------------------------
   1 | 010200000002000000000000000091C340000000000000F03F000000000092C3400000000
00000F03F | Record(CHROM=1, POS=10019, REF=TA, ALT=[T])

My PostgreSQL server is running on a pretty big server: 96 cores, 60Gig of RAM. PostgreSQL and PostGIS versions below:

SELECT version();

PostgreSQL 11.2 on x86_64-pc-linux-gnu, compiled by gcc (Ubuntu 4.8.4-2ubuntu1~
14.04.4) 4.8.4, 64-bit
(1 row)

SELECT postgis_full_version();

POSTGIS="2.5.2 r17328" [EXTENSION] PGSQL="110" GEOS="3.7.1-CAPI-1.11.1 27a5e771
" PROJ="Rel. 6.0.0, March 1st, 2019" GDAL="GDAL 2.4.1, released 2019/03/15" LIBX
ML="2.9.1" LIBJSON="0.11.99" TOPOLOGY RASTER
(1 row)

I researched spatial index generation tuning, and from what I found I tweaked some settings before attempting to create the index:

SET work_mem to '2GB';
SET maintenance_work_mem to '2GB';
SELECT pg_reload_conf();
ANALYZE jem2 (var);
VACUUM;

Then I try to create the index:

CREATE index jem2_spatial_index on JEM2 using GIST (var);

As I monitor top on the server, postgres takes 100% of 1 cpu, but uses nearly no memory....  I find this unusual as I expected the index creation to be pretty memory intensive.  The index generation on the other table representing the linestrings in 3 integer columns shows multiple processes being fired up to build the index and sizable memory used briefly.  Also of note, there is effectively nothing else running on the server, no IO or other resource drains.

Something I did not do was:

alter table jem2 alter column var set storage external;

despite reading it improved query performance, it sounded like it would add time to index creation.  I wanted to factor it out of this slow index creation problem first.


texts: How to adjust spacing and alignment between legend symbol and text in QGIS

texts:
I'm getting some weird and seemingly random spaces between the symbol and text labels in the legend for the Qgis print composer. How can I set the symbols and text to align left/right or divide the symbols into one centered column and the text in another?

Also, in the legend currently, the third entry, "10m interim..." is not displaying the symbol at the proper size. It should be the largest yellow circle.


texts: Topology rule for geometry in ArcGIS Pro

texts:
Is there anyway to write topology rule for the polygon geometry?

I want nice rectangle, without any self-crossing lines. But sometimes I found those. So I want to use topology. I am adding an pic with bad example.

I am using ArcGIS Pro.


texts: Geocoding matches by Street Address or by Zip Code?

texts:
I'm trying to learn more about Geocoding in ArcGIS 10.1 for Desktop.

I've attempted to geocode two ways, both by Street Address and Zip Code.

In the end, the Street Address method yielded more matches, but not by many.

Would anyone be able to give me a quick answer on what type of weakness there is by choosing one method over the the other?

More specifically why Zip Code is the inferior choice.


texts: How can I disable snapping on one feature layer, among many feature layers on the map?

texts:
I'm using ArcGIS Flex API 3.1. I'd like to know if it is possible to choose the layer(s) on which I can snap the graphic I'm moving. Currently, when I activate snapping, my graphics are snapped on all feature layers of the map. I don't want this behaviour.
I want to be able to snap only on one feature class.

Resgards


texts: python-qgis was not found while installing QGIS on Linux Mint

texts:
I am trying to install QGIS on my Linux Mint 20.1 Cinnamon, using this:
https://freegistutorial.com/how-to-install-qgis-on-linux-mint-19/
But when trying to
sudo apt-get install qgis python-qgis qgis-plugin-grass

it says python-qgis wasn't found.
I tried sudo apt-get install qgis python3-qgis qgis-plugin-grass
it showed it depends on following packets :
libgdal20 libqgis-3d3.16.3 libqgis-analysis3.16.3  libqgis-app3.16.3 libqgis-core3.16.3 libqgis-gui3.16.3 python3-qgis qgis-providers qgis-plugin-grass grass740 libgdal20 libqgis-app3.16.3 libqgis-core3.16.3 libqgis-gui3.16.3 libqgisgrass7-3.16.3 .
I tried install them all together with apt-get install it showed
libgdal20 , grass740 ,libgdal20  (no installation candidate exits ) what is that now , I deleted those 3 packets, it showed me another dependecy problems.
here is for  sudo apt install qgis qgis-plugin-grass  :
You have requested an impossible situation or when you have the
Unstable distro use that some required packages still
have not been created or have not yet left Incoming.
The following information may help you resolve the situation:
The following packages have unfulfilled dependencies:
 qgis: Depends on: libgdal20 (> = 1.8.0) but cannot be installed
        Depends on: libqgis-3d3.16.3 should not be installed
        Depends on: libqgis-analysis3.16.3 should not be installed
        Depends on: libqgis-app3.16.3 should not be installed
        Depends on: libqgis-core3.16.3 should not be installed
        Depends on: libqgis-gui3.16.3 should not be installed
        Depends on: python3-qgis (= 1: 3.16.3 + 28bionic) but should not be installed
        Depends on: qgis-providers (= 1: 3.16.3 + 28bionic) but should not be installed
 qgis-plugin-grass: Depends on: qgis-provider-grass (= 1: 3.16.3 + 28bionic) should not be installed
                     Depends on: but grass740 cannot be installed
                     Depends on: libgdal20 (> = 1.8.0) but cannot be installed
                     Depends on: libqgis-app3.16.3 should not be installed
                     Depends on: libqgis-core3.16.3 should not be installed
                     Depends on: libqgis-gui3.16.3 should not be installed
                     Depends on: libqgisgrass7-3.16.3 should not be installed
E: Problems cannot be corrected, you have defective packages

when even try install these dependecies they count on others


texts: OpenLayers, valid zoom levels are offset by 1 level. Why?

texts:
I have an OpenLayers based web app, map object is initialized like this:

p4326 = new OpenLayers.Projection('EPSG:4326');
p3857 = new OpenLayers.Projection('EPSG:3857');

//initialize map object, bind it to map_canvas div
map = new OpenLayers.Map('map_canvas',
    {
    //allOverlays: true,
    projection: p3857,
    displayProjection: p4326,
    maxExtent: new OpenLayers.Bounds(-20037508.3427892,-20037508.3427892,20037508.3427892,20037508.3427892),
    numZoomLevels:20,
    maxResolution:156543.0339,
    units:'m',
    projection: p3857});

I have a pre-rendered map from levels 8-16, however, in order to view zoom level 16 of the map, I have to initialize its layer like this:

var osm = new OpenLayers.Layer.XYZ("OSM Mapnik", "http://localhost/osm_tiles/osm/${z}/${x}/${y}.png", { isBaseLayer: true, sphericalMercator: true, numZoomLevels: 17});
map.addLayer(osm);

Notice, that I set the zoom level 17 (or I can't zoom and view level 16 of the map). Can someone tell me why?


texts: UpdateCursor giving ValueError: need more than 1 value to unpack?

texts:
I am writing a script assigning the ADMINAREA value to the empty LOCALITY fields from TWB_Suburbs to TWB_Property layer. I am getting:

  ValueError: need more than 1 value to unpack

Is anyone able to help out?

import arcpy
# set overwrite true so that if data already exists it overwrite
arcpy.env.OVERWRITE = True
arcpy.env.overwriteOutput = True

# make query for empty LOCALITY values
query = '"LOCALITY" = \' \''

# assign workspace
arcpy.env.workspace = "C:\Users\data"

# make feature layer from property shape file, apply empty locality values query
arcpy.MakeFeatureLayer_management("TWB_Property.shp", "TWB_Property_Layer", query)

# make feature layer for suburbs shape file
arcpy.MakeFeatureLayer_management("TWB_Suburbs.shp", "TWB_Suburbs_Layer")

# Select empty properties in Properties layer WITHIN Suburbs layer
arcpy.SelectLayerByLocation_management ("TWB_Property_Layer", "WITHIN", "TWB_Suburbs_Layer")

#create empty list
suburbList = []

# iterate over admin layer field values, to get it's value
with arcpy.da.SearchCursor("TWB_Suburbs_Layer",["ADMINAREA"]) as rows:
    for row in rows:
        suburbList.append(row[0])
        print (row[0])

#get the number of selected land parcels
numberSuburbs = len(suburbList)

print(row[0])

# Create Update Cursor to update missing LOCALITY fields in property layer
for parcel in range (0, numberSuburbs):
    with arcpy.da.UpdateCursor("TWB_Property_Layer", "LOCALITY") as initialRows:
        for initialRow, rowTwo in initialRows:
                print(initialRow)
                print(rowTwo[0])
                rowTwo[0] = row[0]
                initialRows.updateRow(initialRow)

# Print number of fields updated for each suburb
    print ('Suburb, Number of updated properties' .format(row[0], row[1]))
    print ("{row}" + str(numberSuburbs) +
           .format(row[0], row[1]))

    #report success message
    arcpy.AddMessage("Well done!")


texts: Rasterize(Convert Points to Raster) :UnicodeDecodeError in QGIS3

texts:
I used Rasterize in QGIS3 to convert points to rasters. The source file of the point layer used is as follows:

But clicking Run displays the following error message:
Traceback (most recent call last):
File "D:/QGIS/apps/qgis/./python/plugins\processing\algs\gdal\GdalAlgorithm.py", line 137, in processAlgorithm
GdalUtils.runGdal(commands, feedback)
File "D:/QGIS/apps/qgis/./python/plugins\processing\algs\gdal\GdalUtils.py", line 108, in runGdal
for line in proc.stdout:
File "D:\QGIS\apps\Python37\lib\codecs.py", line 322, in decode
(result, consumed) = self._buffer_decode(data, self.errors, final)
UnicodeDecodeError: 'utf-8' codec can't decode byte 0xb2 in position 12: invalid start byte


texts: Relationship between Android bearing and iOS course?

texts:
I'm a mobile developer working on an Android and iOS app that displays the user's location and various POIs on a map.  I've been comparing the location objects that I get back from Android and iOS and I noticed that they both provide a property that represents the "direction that the device is traveling".  Android refers to this value as a bearing and iOS refers to it as course.

Here's what the Android documentation says:

Bearing is the horizontal direction of travel of this device, and is not related to the device orientation. It is guaranteed to be in the range (0.0, 360.0] if the device has a bearing. If this location does not have a bearing then 0.0 is returned.

And from the iOS documentation:

Course is the direction in which the device is traveling, measured in degrees and relative to due north.  A negative value indicates that the course information is invalid.

Given these definitions, it sounds like they are the same thing but I'm confused as to why they have been given different names.  Is one term (bearing or course) "more correct" than the other or is this just a matter of Google choosing to call it one thing and Apple choosing something else?


texts:
I am trying to  change Qgis splash screen. I edit QGISCUSTOMIZATION3.ini and put this in
[Customization]
splashpath=C:\\GIS\\
...

In folder C:\GIS\ i put png image splash.png, but when I start qgis splash screen not change from default.


texts: Keep nonintersecting polygons during clip

texts:
I am using ArcMap 10. I have two datasets, one of which consists of many small contiguous polygons, and one which has many small non-contiguous polygons. The first dataset is census blocks, and the second is a polygon outline of urban areas generated from the raster datasets NLCD provides. What I want to do is clip each census block such that each census block polygon only contains the urban area specified by the NLCD polygon. Running clip with the census blocks as the input features gives me this. However, there are additional constraints: If the area of the intersection between a census block and the NLCD polygon is under a certain threshold (say, 1000 sq ft), then the entire census block should be kept in the output.

In this example, the green polygon is the NLCD data, and the red outlines are the census blocks. You can see the resulting clip.

The issue comes when there is no intersection between the census block and the NLCD polygon; nothing is kept in the output for that census block.

I need to do this for a large amount of data, probably in arcpy. What strategy would you suggest?


texts: Why are some CAD lines invisible in ArcMap?

texts:
Why am I unable to see certain line features in certain DWG or DXF files in ArcMap?

I've tried changing the symbology, converting the files (with Autodesk TrueView) to various iterations of the DWG standard and even getting a DXF.

In all cases, the lines do not show up in Arcmap at all. I can see the other layers in these DWG files, such as annotation and grid ticks. The most frustrating part is that the missing line is clearly visible in Autodesk.

Caveats: I'm not a trained CAD user. These files are being sent to me for review in GIS software by engineers using (presumably) CAD software.

Conditions: unknown DWG / DXF files. ArcMap 10.2.

Examples: A file that has been giving me this issue.
Screenshot of that file as rendered in ArcMap, with an 80% transparent ortho basemap

Screenshot of that same file in Autodesk Trueview 2014


texts: Why MOLUSCE plugin in QGIS predict the output same as final layer

texts: I have temporal landuse data for the years of 1980 and 1990, spatial variables are Distance to roads, Slope, elevation and Aspect. But the MOLUSCE plugin simulation results with 1 iteration are same as the final layer.


texts:
I was wondering how to avoid the Geoserver transformation that takes place once you request an image via GetMap by using image/geotiff at the end of the url sentence in the browser.

I am experimenting with a geotiff that has 300 x 300 pixels (290 KB) in the server.

If I set the parameters in the GetMap request (url):width=1012&height=1012
The output image will have 1012 x 1012 pixels (1,002 KB). If I request width=300&height=300 - as in the original - the output image will have  300 x 300 pixels and 290 KB.

If I overlap the two images, both have the same bounding rectangle but obviously the pixels don´t match, understood.

But the most intriguing and worrying thing is that if we analyze these 2 output images by using gdalinfo -stats in the command line, you will see that the 1012x1012 image has its parameters totally wrong (Minimum, Maximum, Mean, etc.) but the 300x300 one too.

In fact the original image is a Float32 raster and the 2 output ones are type Byte.

Is there any way of getting the original one from Geoserver, via url, without any transformation at all?


texts: Change countries boundaries color in a shapefile with R

texts:
I am working with R to create a map with a shapefile.
I would like to change the color of the countries boundaries but I don't know how to do that.
Is it possible to change the color?

Here is my code:

world1 = readShapeSpatial("map.shp")
plot(world1)

I have a map with black boundaries and I would like to have grey boundaries.
If I write plot(world1, col="grey") it colors the countries and not the boundaries.


texts: Why is SymbologyControl missing symbols in XP?

texts:
I'm using the the code form SymbolForm.vb in this sample: Change layer rendering using the SymbologyControl It's working great on my developer machine with Windows 7, 64bit.

But when I deploy the solution to an Win XP SP3 32bit the symbols in the SymbologyControl are missing (Text is shown and the symbol can be choosed) when I show the Form in my app.

If I use the designer and add a new SymbologyControl to the MainForm, it shows the same *.serverStyle file correctly.

The problem is caused by the "Enable XP visual styles" in the visual studio project application setting:

When I disable this setting the symbols are shown but not when this is enabled. This is in VB and the code effected by this setting is Application.EnableVisualStyles Does any one know a way round this problem? I don't want to have to disable this on my main project. Actually what's happening is that System.Windows.Forms.Application.EnableVisualStyles() gets runned in the WindowsFormsApplicationBase class.

Update:
Been able to pin this issues down a bit further, this problem only occurs when I show the symbol dialog from events in .Net for example a Button_Click event, when using the events on ESRI components for example OnMouseDown on the TocControl the dialog shows up correctly with the symbols.

I didn't realize it but the problem problem also appears in windows 7, but that issue it's not as big, it only causes the symbols to gets shown in my case two instead of three columns. That's not ok but at least all symbols are shown.

I still haven't found any working work-around, having to remove the call to EnableVisualStyles() is not an option.


texts:
I had a map file for each state (raster file) of India. Now I want to group all map files to show  a single raster or single map file?

How I can do in using QGIS or any platform?.


texts: Get distances within polygons based on a line

texts:
Setup:

Country regions are stored as Polygons in postgres/postgis (Region:id:region_poly)
Cross-country routes are stored as Lines in a different table (Route:id:route_line)

Problem:

Given a Line (route_line), how can I pull the length of the line in each polygon (region_poly) with a single query

Example:

Given Polygon A, B, and C, and a line X that crosses all three of these polygons, I'd like to get something like: {"A": 10km, "B": "160km", "C": 16km}

I think I have to use the ST_Intersection with ST_Length, but I can't quite wrap my head around it. (I'm still new to the gis stuff)


texts: QGIS world map projection with solid Russia (10 degrees over)

texts:
I'm fairly new to QGIS. I have a world map, standard WGS84 with 10 degree graticules (see map).

How do I "shift" the whole map over just 10 degrees to give me a complete world map with a solid (unbroken area) Russia?

See red line on map.

I don't need the size or anything else to change, I just want it to be centered 10 degrees over.

[World][1\


texts: Possible workflows for shaded relief representation as vector

texts:
I am looking for a good workflow on how to represent relief on a map as vector (e.g. for further processing in Adobe Illustrator or Inkscape). My inspiration is the present google maps at zoom level 12, since this looks like it is based on vector graphics:

https://maps.google.de/?ll=47.27853,10.686607&spn=0.179816,0.41851&t=m&z=12

Starting with some free SRTM elevation data, what would be the best way to achieve such a result?


texts: Building a Mosaic Dataset with rasters that are not georeferenced. Whats the workflow?

texts:
I have 68 aerial jpegs that I would like to use in a Mosaic Dataset. They have not been georeferenced. My question is in regards to whether I should georeference them and then rectify and save them in a new folder where I will add later to the Mosaic Dataset, or create a geodatabase, drag all the jpegs into the geodatabase, georeference (rectify?), and then create a Mosaic Dataset inside the same geodatabase.

I noticed that the original file size of the jpegs is around 30 mbs and when rectified creates a 500 mb Tif! Will these sizes be too large to work with when added to a Mosaic Dataset?

I'm using ArcMap 10.1.


texts: I want to change using C# the style of a set of shapefiles once I have uploaded them or while I am uploading them

texts:
I am using C# and the WebRequest class to create a worksapce and then another WebRequest to create a datastore and upload a folder full of shapefiles to a geoserver. However, I would also like to change the style of the uploaded shapefiles to another style that has already been installed on the geoserver. Any ideas?

Lazaros.


texts: Setting new value into field created by PyQGIS?

texts:
I'd like to create a new field named "vctcn93", loop the features of active layer, if string "metro" in feature name, I'd like to add a value string "get" into the field "vctcn93".

But when I check the table, I found a field named "vctcn93" but null value in it.

When I check some features in console, the features got the value already.

These are my codes first time:

 def AddAttri():
    layer = qgis.utils.iface.activeLayer()
    pr = layer.dataProvider()
    feild = QgsField("vctcn93", QVariant.String, "text", 254)
    pr.addAttributes([feild])
    layer.updateFields()

    layer.selectAll()
    ids = layer.selectedFeatureIds()
    features = (f for f in layer.getSelectedFeatures())

    index = a.fieldNameIndex("vctcn93")

    for a in features:
        if "metro_station" in a.attribute("name"):
            a.setAttribute(index, "get")

    layer.updateFeature(a)

And I tried another way to figure it out, but failed.

Field with blank value in the table the same.

Codes second time.

def AddAttri():
    layer = qgis.utils.iface.activeLayer()
    pr = layer.dataProvider()
    feild = QgsField("vctcn93", QVariant.String, "text", 254)
    pr.addAttributes([feild])
    layer.updateFields()

    layer.selectAll()
    ids = layer.selectedFeatureIds()
    features = (f for f in layer.getSelectedFeatures())

    index = a.fieldNameIndex("vctcn93")

    for a in features:
        if "metro_station" in a.attribute("name"):
            fs = a.fields()
            fs.append(feild)
            a.setFields(fs)
            a.setAttribute(index, text)

    layer.updateFeature(a)


texts:
I am trying to use the Table Manager to rename a field in an attribute table.  However, the Vector Toolbar in my install of QGIS 2.14.0-Essen 64-bit Windows shows a toolbar with no tools in it.  The toolbar handle is visible, but not tools are included in it.  See photo below.

Furthermore, when I go into the Customization menu to try and add options to appear in the Vector Toolbar, I see that there are no options to add.  All the other toolbars have a plus next to them that give you the option to add or remove tools.  The Vector Toolbar has no plus.


texts: Create a PostGIS topology with undefined faces?

texts:
I want to have a query that will insert a new edge into an existing topology, and create new nodes accordingly.  I would rather not rebuild the entire topology because I want all the existing node and edge ids to stay just as they are (unless they must be split).  It seems to me that I should follow this workflow:

--SELECT topology.CreateTopology('network_topo', 31467);
--SELECT topology.AddTopoGeometryColumn('network_topo', 'public', 'existing_network', 'topo_geom', 'LINESTRING');
--UPDATE existing_network SET topo_geom = topology.toTopoGeom(geom, 'network_topo', 1, 1.0);
SELECT topology.AddEdge('network_topo', ST_LineMerge(public.new_lines.geom))
FROM public.new_lines;

(The commented out parts are things I have already done)

However the topology.AddEdge command is returning this error:

ERROR:  AddEdge can only be used against topologies with no faces defined

I don't really care about faces, I am only interested in the linear features and nodes.  But I don't know how I can create a topology with undefined faces.

Or maybe I am just on the wrong track entirely.  Is this even possible? Or is it necessary to rebuild the entire topology with a new table of polylines with the new lines added?


texts: Open popup with click on tooltip in Leaflet

texts:
I'm trying to open a popup when clicking on a tooltip of a layer.

I try following code:

function add_layer_rescue_spot() {
    rescue_spot_layer = new L.GeoJSON.AJAX(rescue_spot_path, {
        style: label_style,
        pane: 'label',
        onEachFeature: on_each_rescue_spot_feature,
        pointToLayer: function(feature, latlon) {
            return point_to_label_layer(
                    '',
                    latlon, 'rescue_spot');
        },
    }).addTo(map);

    rescue_spot_layer.on('data:loaded', function() {
        rescue_spot_layer.getLayers().forEach(function(feature) {
            tooltip = feature.getTooltip().getElement();
            tooltip.addEventListener('click', function(e) {
                console.log('test');
                feature.openPopup();
            });
            tooltip.style.pointerEvents = 'auto';
        });
    });
}

function on_each_rescue_spot_feature(feature, layer) {
    rescue_spot_popup(feature, layer);
}

function point_to_label_layer(text, latlon, class_name) {
    var label = String(text);
    return new L.CircleMarker(latlon, {
        radius: 1,
    }).bindTooltip(label, {
        permanent: true,
        direction: 'center',
        className: class_name + '-label',
    }).openTooltip();
}

function rescue_spot_popup(feature, layer) {
    var properties = feature.properties;
    var max_height = parseInt($('#map').height() * 0.95);
    var popup_html = `

Test

                    foo
                    foo

                    bar
                    bar

`;
    layer.bindPopup(popup_html, {maxHeight: max_height});
}

So when I click on the tooltip, 'test' is logged and the map moves a bit like the popup would be opened. But the popup doesn't show up.

When I open the popup manually in the browser by typing:

rescue_spot_layer.getLayers()[0].openPopup()

the popup opens as expected.

Any ideas what I'm doing wrong?


texts: Editing attribute table in ArcGIS Explorer?

texts:
I currently need to edit and add data to attribute tables in Explorer.

Is this even possible?

I am new to Explorer but am proficient in ArcMap. I know how to edit data with ArcMap using the editor feature. The table in explorer came from a geodatabase if that helps.


texts: Create Tiles for a raster(georeferenced) to use in leaflet

texts:
I have raster, with each pixel containing some data.
What i want is to split a large image with each pixel's color is determined by its value, to show on leaflet layer.How to go about this ?

raster 'Size is 19839, 17000'


texts: Eliminate sliver polygons, not the enclaves

texts:
I would like to use the 'Eliminate sliver polygons' in QGIS, but only if they have the same attributes.

To be more specific, I have an administrative shapefile that was created with the SEXTANTE tool 'Union', and which contains cities and villages. Some villages are in cities, some others are half-half, and others are out of the cities. Some of the villages have (very small enclaves) in other villages.
The output of the union tool has a lot of sliver polygons, which I would like to merge.. but if I do that with QGIS 'Eliminate Sliver polygon', selecting polygons per area, it will dissolve my enclaves into the wrong polygon. So I would like to know if there is a way to give a condition for the merging other than biggest polygon/longest boundary.

Is there a way to merge only if this biggest polygon/longest boundary has the same attributes?

My current idea is to 'identify' my enclaves by:

calculating each polygon' area before using the union tool
calculating their area after the union
select polygons that:

a. changed of area

b. are below a certain threshold.

Would there be a faster way to do this?


texts: Convert double to string with Calculator

texts:
I would like to convert double to string using the calculator (or ArcObjects) and get the actual number.  Calculator seems to truncate the number at ~15 decimal places but the actual value may be a lot more precise?  The only way to i can see the actual value is when i hit the "Get Unique Values" button in the Select by Attributes dialog.

Using a random number obtained from the list of unique Values I can write a simple query that correctly retrieves 13 records:

Shape_Length = 9.999983012676239e-005

The value that shows in the Shape_Length attribute is 0.000100000

if I write the following query by simpley copying this value, it does not retreive any records:

Shape_Length = 0.000100000

Taking this a step further, if i create a column and convert the double values from the Shape_Length field into it, though the numbers apear more precise, they are not identical to the value obtained from the Values list.

When I convert 0.000100000 to a string field using the calculator i get 9.99998301267624E-05

When I query Shape_Length = 9.99998301267624E-05 i do get 13 records in this situation but in other cases i may get less records.

in summary:

Shape_Length value: 0.000100000

Converted to string from Calculator: 9.99998301267624E-05

From Get Unique Values list: 9.999983012676239e-005

that's 3 different numbers. I would like to be basically able to see the same number I see in the Unique Values list in my attribute table. How?


texts: Can anybody explain how to compute mean coordinates?


texts:
I am using QGIS 2.6 Brighton.

I am trying to run a Gavrilovic model in the modeller in QGIS.
FYI: Gavrilovic model is a set of Imperical equations for sediment yield analysis.

The equations are as follows:

W = T * h * [pi] * sqrt(z^3)      ..1
T = sqrt ( t/10 + 0.1 )           ..2
z = x * y * ( [phi] + sqrt (J) )  ..3

What I have in QGIS ?

I have raster layers of t, h, x, y, [phi], J.

What do I need to do ?

Calculate z as in equation 3; I have all required variables as rasters.
Calculate T from t as in equation 2.
Calculate W.
Now I need to change the values of t and h rasters. Lets say for example, the existing values of t need to be first changed by -15%, -10%, -5%, +5%, +10%, +15%. And run the equations again to evaluate what will be the changes in W with respect to change in t. Similar for h.
Again I need to evaluate W for -15%, -10%, -5%, +5%, +10%, +15% change in h for each change in t.

  To elaborate point no. 4 and 5.
  Let us assume each raster has 4 cells. And these be the initial values of t and h.

t =  10,  12,  15, 12
h = 110, 100, 140, 90

  Now for the sensitivity test, I need the values of t changed by -15%, -10% , -5%, +5%, +10%, +15%. And similar for h.   6 new t rasters will be created according to the percentage change and the values will be as follows for t.

col A      B      C      D     E     F       G
initial_t -15%   -10%   -5%   +5%    +10%    +15%

   10      8.5     9    9.5   10.5   11      11.5
   12      1.7    1.8   1.9   2.1    2.2     2.3
   15      12.75  13.5  14.25 15.75  16.5    17.25
   12      10.2   10.8  11.4  12.6   13.2    13.8

  Similarly 6 new h rasters will be created, as follows:

col  1       2       3        4        5      6      7
initial_h   -15%    -10%     -5%      +5%    +10%   +15%

110         93.5    99       104.5    115.5  121    126.5
100         85      90          95      105  110    115
140         119     126        133      147  154    161
 90         76.5    81        85.5     94.5   99    103.5

Now, I need to run equations 1 to 3 taking t and h  in following fashion:

for p = col A to col G {
 for q = col 1 to col 7 {
     run equation 3
     run equation 2
     run equation 1
 }
}

  Note:
  I've heard changing values of raster and running this algorithm is easier in graphics modeller.
  If not, I can create those 6 new rasters for t and h myself. I can use R for that. However if it can be done in QGIS (other than using raster calculator) it would be great if you can share the idea.

I need someone to tell me how can I set up a Graphics modeller
and perform  these operations.

  I am not familiar with python but I can find my way writing the code.


texts: Selecting every image of collection using Google Earth Engine?

texts:
I am currently working on google earth engine and I have created an image collection with 15 images in it. What i want to do is to create an image for every single image in my collection. I was able to create an image for the first image of my collection by doing :

var Image1 = ee.Image(Collection.first());

Although i would like to do that for the 15 images. Is there any way to do a loop that reads the entire collection and creates images by iterating. The first one would be Image1, and then Image2, Image3 and so on...


texts: Adding a group layer to another group layer in an Mxd using python

texts:
I would like to take an empty group layer and add it an existing group layer within an MXD.
The empty group layer does not currently exist , i am going to use an empty group layer i have saved within a folder.

i am using a list to build the layers i would like to create and then add.

The whole process consists of having order 1 groups , that contain sub groups.

i was able to do the first step , by creating the intial layers using the code below

Group_1 = ['Scotland','England','Wales','Spain']

for first_Group in Group_1:

        mxd = arcpy.mapping.MapDocument(r"C:\GIS\PROJECT\PYTHON\GROUP_LAYER_TEMPLATE_copy.mxd")
        df = arcpy.mapping.ListDataFrames(mxd)[0]
        groupLayer = arcpy.mapping.Layer(r"C:\GIS\PROJECT\PYTHON\\An_empty_Group_Layer.lyr")
        arcpy.mapping.AddLayer(df, groupLayer, "BOTTOM")
        layers = arcpy.mapping.ListLayers(mxd)
        for layer in layers:
            if layer.isGroupLayer and layer.name == "An_empty_Group_Layer":
                layer.name = first_Group
        mxd.save()

Using the above Code , i was able to add the follwing group layers to the MXD

Scotland

England

Wales

Spain

I Like to use second list i have to create the second level of groups

 Group_2 = [['Scotland','Aberdeen'],['Scotland','Glasgow'],['Scotland','Inverness'],['England','london']['England','bristol']['England','NewCastle'],['Wales',Cardiff],['Wales','Newport']

,['Wales','Swansea'],['Spain','Mardid],['Spain','Valencia'],['Spain','Barcelonia'],

i can use a loop to go through my Second group

That would allow me to add groups , i would need to find a way using the group name to tie into adult group

# An example item in the list ['Scotland','Aberdeen']
# The first element can be used track by to the group above. In this example it would return 'Scotland'. I would like to take this value and then use it to reference the group above.

for item in group_2:
    # x equals the first element of the group
    x = item[0]
#Can use the Second Value from the Item to be used to make the group name
for item in group_2:
    # x Equals Aberdeen
    x = item[1]

Another question related to group layers


texts: Calculating Zonal Stacks on imagery with missing pixels

texts:
I'm trying to calculate zonal statistics on some imagery that has a bunch of pixels that come out as 0s. This is affecting the mean value for the area, or at least seems to be.

I'm using the rasterstats library. Is there some flag to ignore a certain value?


texts: OpenLayers plugin OverflowError: cannot convert float infinity to integer

texts:
OpenLayers plugin in QGIS 2.4 is not working properly when not loaded in first place.

If I add an OSM map with the plugin immediately when opening QGIS, all goes well. But  when I call the plugin after having added some other layer before, I have these trace-backs:

An error has occurred while executing Python code:

Traceback (most recent call last):
File "/home/username/.qgis2/python/plugins/openlayers_plugin/openlayers_layer.py", line 134, in pageLoaded
self.setup_map()
File "/home/username/.qgis2/python/plugins/openlayers_plugin/openlayers_layer.py", line 184, in setup_map
olSize = QSize(int(olWidth), int(olHeight))
OverflowError: cannot convert float infinity to integer


texts: How can I compute raster pixel width and height given raster bounds, row count, and column count?

texts:
Restated, how do I convert arc degrees to meters? For example, I have an elevation data raster that has the following metadata:

WEST LONGITUDE=64.97798789° E
NORTH LATITUDE=33.02003415° N
EAST LONGITUDE=66.03338707° E
SOUTH LATITUDE=31.98030163° N
PROJ_DESC=Geographic (Latitude/Longitude) / WGS84 / arc degrees
PROJ_DATUM=WGS84
PROJ_UNITS=arc degrees
EPSG_CODE=4326
NUM COLUMNS=9001
NUM ROWS=9001
PIXEL WIDTH=0.0001111 arc degrees
PIXEL HEIGHT=0.0001111 arc degrees

I can compute by hand the pixel width in arc degrees as follows:

(EAST LONGITUDE - WEST LONGITUDE) / NUM COLUMNS

Similarly, I can compute by hand the pixel height in arc degrees as follows:

(NORTH LATITUDE - SOUTH LATITUDE) / NUM ROWS

My question is how to compute the pixel width and height in meters.


texts:
I have a lot of rasters. Here is the typical gdalinfo output:

Size is 14058, 9940
Coordinate System is `'
Origin = (1521634.513864640844986,403292.969188707182184)
Pixel Size = (0.169386038323864,-0.169386038265657)
Metadata:
  TIFFTAG_DATETIME=2019:11:18 08:53:19
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=Adobe Photoshop CS6 (Windows)
  TIFFTAG_XRESOLUTION=72
  TIFFTAG_YRESOLUTION=72
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( 1521634.514,  403292.969)
Lower Left  ( 1521634.514,  401609.272)
Upper Right ( 1524015.743,  403292.969)
Lower Right ( 1524015.743,  401609.272)
Center      ( 1522825.128,  402451.121)
Band 1 Block=14058x1 Type=Byte, ColorInterp=Red
Band 2 Block=14058x1 Type=Byte, ColorInterp=Green
Band 3 Block=14058x1 Type=Byte, ColorInterp=Blue

Each file is around 400 mb, which is too much for me. How can I decrease the file size of each raster?


texts: arcpy.Delete_management not deleting a geodatabase folder

texts:
This script runs fine the first time but fails when run the second time. The issue seems to be that the KMLToLayer_conversion statement creates a file geodatabase (no surprise) that subsequently cannot be deleted even when the feature class is removed from the map, the layer file is deleted, and the geodatabase contents are deleted. I'd like clean up after myself when this script is done with all traces gone other then the new feature classes within MasterGDB. The problem is that this script can only be run once unless you exit ArcMap, manually delete the folder in Windows, then re-start ArcMap. Running the individual commands within the Python window inevitably shows a "" but the directory for the geodatabase remains. So what's this newbie missing here?
(In this test/debug script, there is only a single KML file "C:\Temp\KKKLLL.KML"  -- which converts just fine).

import arcpy, os
# Name: BatchKML_to_GDB.py
# Source: AS16818.ZIP from acripts.esri.com

import arcpy, os

# Set local variables and location for the consolidated file geodatabase
KMLDir = "C:\TEMP\KML3"
outLocation = "C:\\Temp\\MuleDeer"
MasterGDB = 'AllKLM5.gdb'
MasterGDBLocation = os.path.join(outLocation, MasterGDB)

# Create the master FileGeodatabase as needed
if not (arcpy.Exists(MasterGDBLocation)):
    print MasterGDBLocation + " doesn't exist; creating it now"
    arcpy.CreateFileGDB_management(outLocation, MasterGDB)

# Convert all KMZ and KML files found in the current workspace
# Set workspace (where all the KMLs are)
arcpy.env.workspace=KMLDir
for kmz in arcpy.ListFiles('*.KM*'):

  print "CONVERTING: " + os.path.join(arcpy.env.workspace,kmz)
  kmz2 = os.path.join(arcpy.env.workspace,kmz)
  arcpy.KMLToLayer_conversion(kmz2, outLocation)
  print "Done"

# Change the workspace to fGDB location
arcpy.env.workspace = outLocation

# Loop through all the FileGeodatabases within the workspace
wks = arcpy.ListWorkspaces('*', 'FileGDB')
# Drop Master GDB from the array/list
wks.remove(MasterGDBLocation)

for fgdb in wks:
  # Change the workspace to the current FileGeodatabase
  arcpy.env.workspace = fgdb

  featureClasses = arcpy.ListFeatureClasses('*', '', 'Palacemarks')
  for fc in featureClasses:
    fcCopy = fgdb + os.sep + 'Placemarks' + os.sep + fc
    arcpy.FeatureClassToFeatureClass_conversion(fcCopy, MasterGDBLocation, fgdb[fgdb.rfind(os.sep)+1:-4])

arcpy.Delete_management(fcCopy)
arcpy.Delete_management("C:\\Temp\\Muledeer\\KKKLLL.lyr")
arcpy.Delete_management(fgdb)


texts:
I want to track and visualize the progress of a bus, as follows:

I have the coordinates of each station
as well as the coordinates of the bus, updated every 10 seconds

The information is stored in Postgres, as float tuples, without any awareness of street topology, or other GIS-related metadata. At the moment I have a naive implementation that calculates the "as the crow flies" distance between 2 points (the bus and the next expected station). This allows me to build a rough progress-bar, but due to the way the streets are laid out in certain parts of the city, sometimes these estimations are... very wrong.

My understanding is that the proper way to deal with this is to rely on "map-matching", which will do a "snap to road" transformation of the coordinates of the bus. This, in turn, will allow me to see what distance was covered, as if the entire route is a straight line:

However, I do not have the GIS experience that would give me the confidence that (1) I am on the right track and that (2) I chose the path of minimum resistance.

What are the best practices in addressing this problem?

Update: the current software is a pure Python program that performs these calculations. No GIS tools are involved at all, though I suspect that the problem at hand is probably a typical scenario for them. I would prefer a solution that uses free/libre open source software, which can be integrated into existing systems, a Python API would be a plus.


texts: Topological relationship between two polygon feature layers in VBA and ArcObjects?

texts:
I have two polygon feature layers: one representing land parcels and one representing roads.

Some land parcels (polygons) coincide with a part of a road (another polygon feature). How can I check which land parcels coincide in any of their boundaries with a road?

Of course help will be great if a code in VBA and ArcObjects is already available.

I note also that I have considered the topology rules defined in http://webhelp.esri.com/arcgisdesktop/9.3/index.cfm?TopicName=Topology_rules and also and the ITopologyRule interface but I could not be sure which rule in both cases to use for solving my problem.


texts: Problem using Semi-automatic classification plugin on QGIS for downloading Sentinel-2 images

texts: I´m trying to use the Semi-automatic classification plugin on QGIS 2.14.3 for downloading Sentinel-2 images. No secret following this video tutorial https://www.youtube.com/watch?v=HoHSlKfSiCE&feature=youtu.be. But I don´t know why, does not work. Error [50]: Internet error.


texts: From several polygons of one layer to one polygon

texts: I have imported one CSV file where I have zones of the city. Is it possible to get one big polygon of the city?


texts:
I ran a script that I got from exporting a model that I built in Model Builder in ArcCatalog 10.3.1 and made edits to it (using the generous help from people on here). I essentially put in for my buffer, feature class to feature class, and spatial join processes to print if they already exist and then the code to delete them.

After running the script, I received this text in blue (I didn't play around with the color defaults on Python IDLE):

F:\Projects\Working\ModelProject\ModelData.gdb\Hurricanes_in_Counties_Feature_Class already exists!!
F:\Projects\Working\ModelProject\ModelData.gdb\Hurricanes_Buffer already exists!!
F:\Projects\Working\ModelProject\ModelData.gdb\Hurricanes_Counties_SpaJoi already exists!!

I was told that I need to put in the output files. My question is how do I even begin with an output file?

Here is my code for the buffer process:

import arcpy
Hurricanes_Buffer = "F:\\Projects\\Working\\ModelProject\\ModelData.gdb\\Hurricanes_Buffer"
if arcpy.Exists(Hurricanes_Buffer):
    print "{} already exists!!".format(Hurricanes_Buffer)
    arcpy.Delete_management(Hurricanes_Buffer)

arcpy.Buffer_analysis(Hurricanes_in_Counties, Hurricanes_Buffer, v15_Miles, "FULL", "ROUND", "NONE", "", "PLANAR")


texts: Georeference a image representing the world with unknown projection system

texts:
I have a set of images like this one attached below:

I would like to georeference the image/map in order to overlay it with Google Satellite images and then produce some polygons in a GIS.
Unfortunately, when I try this operation (through QGIS' plugin 'Georeferencer'), I obtain maps with strong distorsion at borders. For example, it is impossible to overlay the countries of Oceania (which is the main focus of my study).

I think that depends (but I'm not sure) from projections systems, which affect the borders of a map.

Any idea to solve this trouble?


texts:
I'm trying to export a portion of Natural Earth's world 10M hillshade in ArcMap but above a certain latitude, the resolution unexpectedly changes and is really poor.
This is a screenshot of the area in Arc. The resolution over the entire area is fine.

When exported, there is a very distinct line where the resolution abruptly changes.
Hillshade image is imported into data frame (projection is Canada Lambert Conformal Conic). It's being exported from Map View as follows: File--> Export map - exporting as Tiff. I've tried 300dpi to 600 dpi resolutions and the issue remains the same. I've also tried exporting in the same fashion as .jpg. The same thing happens.

Does anyone have any idea what might be causing this to occur? Any wisdom would be greatly appreciated.


texts:
ESRI offers a ParksFinder map that's part of their LGDM (Local Government Data Model) that they have begun offering freely. This was on an external-facting server running 10.1 ArcGIS Server, connected to a 10.1 Enterprise GeoDB, using a 10.2.2 Desktop with the 10.2 version of the ParksFinder. We couldn't get the 'Search by Activity' to work. At all. The other searches gave us issues as well but we discovered little adjustments we could make to the Database, or through the Desktop MXD file, etc. that cleared up all of the searches EXCEPT for the one by Activity.

We wanted to showcase this functionality to our Parks department despite these issues so we setup the same app on an internal-facing server. This server was running 10.2.2 ArcGIS Server, 10.2.2 ArcGIS Desktop, but connected to the same 10.1 Enterprise GeoDB. We have to hack around and change up some minor edits but we got it working.

Chalking the issues up to a difference in Versions of our software and DB we waited for our Quarterly Release when we could upgrade our servers. Which we did yesterday.

So now it was time to implement the ParksFinder app on the newly updated servers. We ported the app over exactly the same. The config.js and proxy.config files were updated with the correct url's to services and everything seemed to work... EXCEPT THE SEARCH BY ACTIVITIES!

We spent 4 hours, adjusting permissions, rebuilding application pools, deleting and re-creating, everything... but nothing worked. I had gotten to work at 6 that morning and by 6:41 that night I was done. My co-worker (a much more advanced GIS professional) and our GIS DBA couldn't let it go and kept trying out as many things as they could.

The thing that fixed it is what has me stumped.

In the locator.js file of an untouched ParksFinder app there is a variable called activityQueryString that creates the URL Query Builder for the searches and around line 1560-ish it get's this little addition.

activityQueryString += "1=1";

I would never have thought to do what they did... by changing that line to the following:

activityQueryString += " 1=1 ";

Adding spaces to the query string fixed the issue.

My questions is, why?

What does that add to the search functionality that would turn it from non-working to working?

I'm confused. I've got a LOT of experience with web dev, mainly with CMS systems and basic websites, and I've dealt with query strings being difficult before but nothing like this.

Just trying to make sense of a crazy situation.


texts: How to individually label points in Openlayers linestring?

texts:
I want to display some text below each point on a linestring with Openlayers.  The text would include, for example, the coordinates of the point, plus some other information.

I have managed to add a style to these points, but so far only the same style to all.

https://jsfiddle.net/wainwrightalex/6qbpqmme/8/

I considered creating a new feature for each point, but I would like to avoid this if possible, so I can still use ol.interaction.Modify to edit the linestring.

Is it possible to somehow style these points individually?  Or add individual text in some other way?


texts: Creating custom scale measurement for QGIS project

texts:
I started to make a map based on a raster (a scanned paper map). The map obviously has not geographic coordinates or a projection.
The map comes with its own scale, and I would like to somehow make my QGIS project use the same scale.

As you can see, the map uses a scale of 1 inch = 5 linear arpents (the arpent is a French area unit which also has a linear counterpart). 1 linear arpent = 58.471308 meters.
If I use the measurement tool, it says that the inch measures 0.016 meters. How can I set QGIS to accept that from now on, this value shall be 5 arpents ? I would want my scale box in the printed form of my map to be accurate.


texts: Centering polygon on specific coordinates using QGIS

texts:
I need to place a set of polygons onto a set of points, so that each polygon is centered on its corresponding point.
How do I do that?
In MapInfo I can just manually edit the coordinates of the center, but there does not seem to be an equivalent in QGIS.


texts: Validating schema of my XSD wth GML with GeoTools

texts:
I would like to control if my GML file is valid with my actual XSD using GeoTools.
I found something here: GML2 Parsing TestSchema Example. But I don't know if it's a good way to do what I want.
Configuration configuration = new org.geotools.gml3.GMLConfiguration();
Parser parser = new org.geotools.xsd.Parser(configuration);

//the xml instance document above
ClassLoader classloader = Thread.currentThread().getContextClassLoader();
InputStream xsd = classloader.getResourceAsStream("gml/config.xsd");

int i =0;
//parse
SimpleFeatureCollection fc = (SimpleFeatureCollection) parser.parse(xsd);
fc.accepts( new AbstractFeatureVisitor(){
    public void visit( Feature feature ) {
        SimpleFeature f = (Feature) i.next();

        Point point = (Point) f.getDefaultGeometry();
        String name = (String) f.getAttribute( "name" );
    }
}, new NullProgressListener() );

In my code I don't know where I can specify my GML file for validation.


texts:
My code runs like this:

arcpy.gp.RasterCalculator_sa("""Int(RoundUp("time_min"))""",path+"/Time_int")
arcpy.ExcelToTable_conversion("D:/Flash Flood/Python/Reclassify.xlsx",path+"/reclassify","5")
arcpy.gp.ReclassByTable_sa("Time_int","reclassify","From_","To","New",path+"/Rectime","DATA")
Classify = arcpy.GetParameterAsText(12)
arcpy.ExcelToTable_conversion("D:/Flash Flood/Python/Reclassify.xlsx",path+"/reclasviz",Classify)#error
arcpy.gp.ReclassByTable_sa("Time_int","reclasviz","From_","To","New",path+"/Viz","DATA")

For some reason the second reclassify doesn't work. Why is it so? The input itself runs into a problem. This is only when I run the whole program. If I just separately run the last 3 lines providing the layers, it works fine. The Excel file has classification limits.


texts: Embedding interactive map maker in external website

texts:
I'm currently developing a new website that will allow users to submit project proposals via an online form for assessment and feedback. The system is currently working well, but I'd like to extend it to allow people to submit spatial information as part of their project submissions (it's for global landscape projects). Within the current form, I'd like users to be able to create maps either from fresh or via upload of existing spatial data that they can then save and submit.

I have access to ArcGIS Online and ideally would like to embed the map creation functionality within my form on my external site that is then saved and submitted, however I'm not sure if this is possible?

I've read through the API documentation but couldn't spot anything that describes embedding in external sites. I'm also open to using alternatives such as embedding Google Map Maker (if possible?), but I'm guessing this will run into problems where each form user will require a Google account and be pre-authenticated somehow before submitting the form (a big disadvantage).

Does anyone have any ideas on how I could achieve such a task?


texts: Extracting stratigraphic profiles from multipatch features representing continuous geological formation in ArcGIS Pro?

texts:
I have some multipatch undergroud layers representing each a specific geological formation unit. These units are vertically continuous and I'd like to "clip" them with an orthogonal plane to extract a vertical stratigraphic profile.
In order to do that, my simple idea was to:

draw one or more 2d lines representing the profile line seen from above.
extrude this by negative heights so that it goes through all the units.
convert the extruded line to multipatch.
merge the geological unit multipatches together.
use the multipatch from point 3. and 4. as input for Intersect 3D.

So far, I was able to reach the final point, but I receive a warning message reading "WARNING 050084: No Intersect result generated", and no feature is returned. I guess the problem may arise from the multipatch in point 3. being "width-less" (flat) because derived from an extruded line. I cannot even close it, if this could be a problem. Below is an image to clarify what my data looks like (only two geological units are shown in this case).

Does anyone have any other idea on how to achieve this, or how to solve the warning?


texts: Splitting line GeoDataFrame basing on given lengths

texts:
I am trying to split line GeoDataFrame into sections. I have lengths on which splits should be made stored in column. My idea was to generate points along the line using shapely.ops.interpolate and then split lines with shapely.ops.split and collect all the results in new GeoDataFrame.
Unfortunatelly, shapely.ops.split returns given line as if generated points would not intersect with it.
Here is my code:
lines = []
for split_distances, line in zip(suln1234['split_distances'], suln1234['geometry']):
    points = []
    if split_distances:
        for split_distance in split_distances:
            points.append(line.interpolate(split_distance))
        splitted = shapely.ops.split(line, shapely.ops.snap(shapely.geometry.MultiPoint(points), line, 1))
        for i in splitted:
            lines.append(i)
    else:
        lines.append(line)

l_split = gpd.GeoDataFrame({'geometry': lines}, crs=suln1234.crs)

As you can see, I also tried to use shapely.ops.snap, but even with crazy big tolerance, it only works for part of lines.
What is the reason for such behaviour? I would imagine that if points are derived directly from line, they should allign nicely with its geometry. How to fix it?


texts: QGIS use multiple rasterfiles: virtual raster with Python

texts:
I'm new to QGIS (using version 3.4). I'd like to make our raster files ready to use. There are 5k rasters in one folder with several subfolders. What's the best way to display them all in QGIS?

I thought that best would be to use a virtual raster.
Because of the many subfolders and because I need to repeat the process of creating the virtual raster, I wanted to write a Python script.
That's how far I've got trying to add only one file to the vrt:

from osgeo import gdal

input_File = r"D:\Folder\Raster.ecw"
output_File = r'C:\Folder\VRTFile.vrt'
gdal.BuildVRT (input_File, output_File)

This doesn't generate a traceback but also no file.
What is wrong?

I also tried this:

import subprocess
subprocess.Popen(['gdalbuildvrt','C:\Folder\VRTFile.vrt','D:\Folder\Raster.ecw'])

with the same result.

(Also, how can I import all files from the folder? Can I use similar expressions as with arcpy?)


texts: How to make sure GDAL Java binding doesn't crash my whole application?

texts: I'm using GDAL Java binding  to convert some files (ogr2ogr), such as shapefile to GeoJSON. The problem is when something goes wrong it crashes the whole application and I have to start it up again. Is there a way to prevent this?


texts: Exporting table of daily temperature (ERROR: Internal Error)

texts:
I am trying to download a table of daily temperature for the year of 2001 and all Brazilian municipalities in Google Earth Engine Code (GEE). After mapping a Reduce Regions function and flatten to convert in a table, I hit Export to Drive, and it runs for a while (it's a lot of information) but I end up with an internal error. I've tried lots of different code but always ending with this error. I don't think it is a problem in my code. Can anyone help me?
// Script for daily temperature and precipitation extraction for all Brazilian municipalities
// 30/06/2020

// Importing Brazilian municipalities feature collection
var munBrazil = ee.FeatureCollection("users/rodrigoccuri/BRMUE250GC_SIR");

// Importing MODIS Terra Land Surface Temperature daily global 2001

var modisTerra = ee.ImageCollection("MODIS/006/MOD11A1").select(['LST_Day_1km']).filterDate('2001-01-01', '2002-01-01');

// Creating function for Image Collection reduction

var reduceTemp = function(imagemodis){
  return imagemodis.reduceRegions(munBrazil, ee.Reducer.mean(), 1000);
};

// Maping the function to all image collections

var munModis = modisTerra.map(reduceTemp);

// Flattening the collection of collections into a table

var tableMunModis = munModis.flatten();

// Printing the fisrt element to check it out

print(tableMunModis.first());

// Exporting to csv file

Export.table.toDrive(tableMunModis, 'temp2001', 'GEE', 'temp2001');


texts: How to rotate text in ArcGIS by any chosen number of degrees?

texts:
I want to add a textline to my map and rotate it by about 160 degrees, but when right-clicking the element and selecting rotate/flip there is only the option to rotate by 90 degrees (100 US degrees). Is there a (non-scripted) way to do this, or perhaps a simple workaround?

The short text shall be put near one of the edges of a polygon, but does not belong to a specific feature that I want to label. I suppose I could make a new feature class and draw a line along the edge of the polygon, and then label it, but I am wondering if there is a better way to do it.


texts: Performing in-line variable substitution for path in ArcGIS ModelBuilder?

texts:
I am having difficulties in finding a way to use the input path and set is as the output path for a simple Copy Features tool. I use the substitution by the workspace but it outputs the results into the default ArcGIS folder.

I am not sure how to set the workspace to hold the path of the input feature class. As the input feature class is a parameter, I cannot use a hardcoded path. Please see below:


texts: How to know the percentage of features that represent a classification in a vector layer

texts:
I have a classified vector grid layer  with five classes

And I want to know what percentage each class represents of the whole vector layer.

I'm currently selecting the cells from each class multiplying by 100 and then dividing by the total number of cells in the grid, is there a faster way of calculating the percentages?


texts: Drawing Map on two or more sheets without Overlays?

texts: I am producing a map for the forest management department. The map is about a forest station and has a strict scale of 1:10000 on AO paper size but my "station" is too large and cannot fit into one sheet/paper. If I produce it on two sheets there definitely an overlap 'where features are repeated in both sheets'.


texts: How to calculate the percentage of the area that a polygon covers a cell in a raster in QGIS

texts:
I'm having trouble with using QGIS.

I have a GRIB file containing meteorological data, this is loaded in QGIS as raster data if I'm not mistaken.

I also have a shapefile containing a polygon, describing regions within the area covered in the GRIB file.

What I need to know, is how many % of a raster-cell is covered by a region within the polygon. I'm using QGIS.

What I have done so far, is I rasterized the polygon, and used zonal statistics to calculate the SUM and COUNT on the shapefile, and then used the field calculator to calculate the percentage covered. This does not seem to provide the result I hoped for.

Can anyone push me in the right direction?


texts: Broken Data Connections between Esri Maps for Office and ArcGIS Online after update?

texts:
We have a workflow in place where a spreadsheet is updated daily with new location data, synced to ArcGIS Online (AGOL), and then viewed in an AGOL web map and web app for informational purposes. Every time AGOL pushes out an update, the data connection between that spreadsheet and AGOL is severed and the app becomes out of date, even though we keep syncing changes from Excel and the layer says that it's been udpated.

Esri has no solution since they can't recreate the issue because the issue only happens in the gray area between updates, even though they have a record that this happens to us literally every time there's an update.

Is there anything I can do to fix this, or decrease the frequency of failure on my end?


texts: Change Circle Size According to Value using R

texts:
How do I change circle size in a population map according to the population density of the location?

library(maptools)
library(RColorBrewer)
library(classInt)
library(maps)

## load the file ##
pop <- read.table("[pop.txt][1]", header=TRUE, row.names=1, sep=",")

## select color palette and the number colors to represent on the map ##
colors <- brewer.pal(9, "YlOrRd")

#set breaks for the 9 colors
brks<<-classIntervals(pop$POP10, n=9, style="quantile")
brks<- brks$brks

map(database='state',regions=c("California"),col="black", lwd=.3, resolution=.1)
points(pop$LONG, pop$LAT, pch=21, bg=colors, cex=0.7, lwd=.4)

This way I obtain a map like this.

I would like to create a map that shows very small dots for little population density (lighter points in the example) and bigger dots as the population density increases.


texts: How to use gdal_rasterize with postgis vector?

texts:
I'm trying to convert a multipolygon vector layer to a raster using the GDAL tool gdal_rasterize. The command I am using is:

gdal_rasterize -b 1 -burn -9999 PG:"dbname=db password=password" \
    auckland auckland.png

However, it seems my specification of the source for the postgis layer is incorrect as I get an

ERROR 4: 'auckland.png' does not exist in the file system,
and is not recognised as a supported dataset name.`

So how do I successfully connect to postgres and export a vector as a tif? Is it possible?

Numerous other attempts have included:

gdal_rasterize -b 1 -i -burn -9999 PG:"dbname=db password=" -l \
    auckland auckland.png

gdal_rasterize -b 1 -i -burn -9999 -l auckland PG:"dbname=db \
    password=password" auckland.png

And this one from Roger's Rants:

gdal_rasterize -b 1 -i -burn -9999 PG:'host=localhost dbname='unep'' \
    -sql "SELECT ST_buffer(the_geom, .12) from coastlines where name='Poland'" \
    Poland_elev.tif

When I use ogr2ogr to export the layer out of postgres as a shapefile, and then gdal_rasterize to convert to a .tif file it all works smoothly.

Using GDAL 1.11.2 version on Ubuntu 15.10.

Still no love with either:

gdal_rasterize -b 1 -burn -9999 PG:"dbname=db password=pwd" \
    -l auckland auckland.tif

or:

gdal_rasterize -b 1 -burn -9999 -l auckland PG:"dbname=db \
    password=pwd" auckland.tif


texts: Create a Menu Containing a List of Layers That Can be Added to an MXD Using arcpy.mapping

texts:
Is it possible to create a list of layers (residing in something like a menu) that can be selected and added to an mxd using arcpy.mapping.  I've been using the bit of code below (there's someone I need to thank for this as I nabbed it from this site) to add a single layer. The code is within a Scripting Tool and is added to a toolbar as a Geoprocessing Tool.

import arcpy
import arcpy.mapping

#Data Location
parcels = r'P:\Mapping\Data\Land.gdb\Parcels\Parcels'

#Get the map document
mxd = arcpy.mapping.MapDocument("CURRENT")

#Get the data frame
df = arcpy.mapping.ListDataFrames(mxd, "*")[0]

#Create new layer
parcelLayer = arcpy.mapping.Layer(parcels)

#Add the layer to the map at the bottom of the TOD in data frame 0
arcpy.mapping.AddLayer(df, parcelLayer, "TOP")

#Refresh things
arcpy.RefreshActiveView()
arcpy.RefreshTOC()
del mxd, df, parcelLayer

I want to avoid creating a script for each dataset that I want to add programatically to an mxd, but I also want to stay within Python/arcpy as I have no business trying to do this within ArcObjects. Any help is appreciated.


texts: QgsfeatureRequest - Store result indexes

texts:
I recently posted a question. My original goal was to find a way to determine on which side of a line points are located. @xulnik answered this question with two codes. One returns a binary [0;1] column indicating on which side points are located. However as the number of points for which I want to know whether they are on the right or the left of a given line is huge (about 1,490,000) @xulnik proposed an alternative based on enumerate list for points features point indexes.

In his answer @xulnik proposes that I use indexes lists to produce a memory layer with QgsFeatureRequest. I have tried things based on @gsherman 's answer to this post but did not succeed. Basically I am trying to find a way to either export the result of the following code to a csv file or to create a layer that contains the result of the code.

from PyQt4.QtCore import QVariant

registry = QgsMapLayerRegistry.instance()

points = registry. mapLayersByName('pixels_id')
line = registry. mapLayersByName('line')

points_feats  = [ feat for feat in points[0].getFeatures() ]
line_feat = line[0].getFeatures().next()

idx = points[0].fieldNameIndex('position')

ids_0 = []
ids_1 = []

for i, feat in enumerate(points_feats):
    pt = feat.geometry().asPoint()
    point_in_line = line_feat.geometry().closestSegmentWithContext(pt)[1]
    az = pt.azimuth(point_in_line)

    if az > 0:
        pos = 0
        ids_0.append(i)

    if az < 0:
        pos = 1
        ids_1.append(i)


texts: OL5 functions in JavaScript not working with NodeJS

texts:
So I made a weathermap with GeoServer and OpenLayers5.

Works fine and all, but NodeJS keeps giving me a hard time. Does anybody know how this really works?

Because I need input from the end-user, but NodeJS doesn't understand that. I need to be able to put in a date and update my WMS params so that I can show different types of weather (call different kind of layer groups from GeoServer). What the OL5 setup tells me to do is this: https://openlayers.org/en/latest/doc/tutorials/bundle.html. (NodeJS with NPM and Parcel/Make an index.html and index.js/Create a bundle/I can develop using NPM start/To create a production bundle I only have to 'npm run build' and then copy those dist files to my production server.) Seems easy enough. I did all the steps, and indeed I can now show a map with my GeoServer layers.

But what I CANNOT do is make a function in my JavaScript and use that for my user-input. Whenever I put a function in my index.js and want to use it in my HTML it tells me that it is not defined.

I asked a sort similar question on Stack Overflow and they directed me to Express. And I spend my whole morning watching tutorials on Express, but I do think it should not have to be this difficult. Especially since OL5 setup guide does not talk about any of that. Who can help me?


texts: Include gridlines with attribute table shown on a map layout?

texts: I am including the related attribute table to the spatial data shown in a map; see attached image. I have the attribute table shown with the map in layout view, but I would like to include gridlines to separate the individual lines of the attribute table. The table is added utilizing the "add table to layout" option available in the attribute table dropdown menu. Although this shows the gridlines in layout view, the gridlines do not export with the pdf option.  How can I include table gridlines in the final PDF map?


texts: Differences between Operations Dashboard for ArcGIS and Power BI?

texts:
We are looking for a possibilty to show some live statistics and calculation over geodata (the developement on some retail markets ) and want to see what is suitable. I have tried Power BI and have seen there is a possibilty to show maps and statistics at the same time.

I have read About Operations Dashboard for ArcGIS too but unfortunately do not have direct access to that to be able to test it.

Can you guide me about the necessity of Operations Dashboard for ArcGIS and its advantages in comparison to Power BI?


texts: Smooth Data Values Across Geographies with Zero Population

texts:
Looking to fill in gaps through generalizing a population count across U.S. counties. I.e., if a county has 0 population, but is surrounded by 5 non-zero counties, then change the 0-population county to an estimate based on the neighboring non-0 counties. This is all vector polygon data, no points, no raster, etc.

I do not have the Spatial Analyst extension in ArcMap or Pro, so I have been trying the R-ArcGIS bridge extension to run a Empirical Bayes smoothing on the data. The issue is that the calculation fails when it encounters non-positive population values.

Another thought was to use the Polygon Neighbors tool which will spit out counts for all adjacent geographies, and maybe from there I could assign the average of the neighbors. Does anyone have a more precise way? Even if it's a Pythonic?


texts: Rectification of a lot of images using ERDAS IMAGINE

texts:
I want to rectify many images of same area but different dates at a time in ERDAS Imagine
1 I did this using Batch processing but I have no out put. Second I did reproject many images using Batch process. It ran but the projection did not change. What is the error?

If it could not be happen in ERDAS IMAGINE then what is the purpose of this Batch in ERDAS Imagine Rectification process


texts: using mod() with IF THEN field calculation

texts:
I am trying to use the field calculator to assign a "odd" or "even" value to a field of small integers.  I usually do this with by selecting where:

MOD([field],2) = 0

and then just calculating it as "even", then reversing the selection and setting it to "odd"

I was hoping to do something like this in the field calculator:

DIM oddeven
IF mod([Field],2) = 0
   oddeven = "even"
else
   oddeven = "odd"
end IF

This does not work and I assume there is an alternate use of the mod() in this setting, or perhaps you cannot use an operator in the IF/THEN condition statement.


texts: Are feature ids in QGIS and ArcMap stored differently?

texts:
I have a shapefile that I use as a layer in ArcMap. Opening the attribute table, shows a fid starting with 0 for the features. If I attach the same shapefile in QGIS and open the attribute table, the features increment from 1. Copying and pasting the table into Excel gives me an FID field for Arc but not for QGIS. However, Copying the layer from QGIS into SQL Server gives me a table with a key field "qgs_fid" starting with 1.

Does this mean that I cannot treat the fid as a key field that is associated with the record for joins and relates, especially when combining with external data? Does the fid or qgs_fid change if records are inserted or deleted?


texts: Adding tile size and compression type for ArcPy code?

texts:
I have to add tilesize and compression type to the following code

import arcpy

arcpy.CopyRaster_management("c:/data/background.tif",
                           "c:/fdgb/CpRaster.gdb/fgdbRD","DEFAULTS",
                           "0","9","","","8_BIT_UNSIGNED")

How do I add these two lines to above code?

arcpy.env.compression ="NONE 75"
arcpy.env.tileSize ="0.1 0.1"


texts:
I know this question has been asked already in this forum, but I couldn't find an easy (for me, a rookie) solution for my problem.
I have an RGB 3 band file which, in some pixels, has i.e.

R no data
G 147
B 147

This creates problems in the visualization of the raster. I'd like to know how to change those values to 0, enabling the colours to be  correctly displayed again.
I tried to use an algorithm in the processing tab, but it says SAGA algorithms are not supported by multiband layers. Also I've seen many using GRASS tools, but I really have no clue on how to use GRASS, so I would need a step by step explanation. Consider I am still a newbie with this software so, if possible, I'd like use the easiest method.


texts: Using "raster.transform" function of "rasterio" in Python

texts:
When I execute the "raster.transform" function in Python it returns me a matrix.

What I can understand is that it tells me my data has a 10 m pixel resolution,
It also tells me that my left bounding box is at 728749.5 and my top bounding box is at 4373200.5. What I don't understand is the rest of the elements like "-10", "1", and "0s". I also noticed that it can be a 2x3 matrix. Why is this dimension different?

I would like to know what does this matrix conveys and how can I use this information, especially when I have this info from two different rasters.


texts: ArcGIS JavaScript API set size of CSV Marker by attribute

texts:
Using ArcGIS Javascript API 3.29 is it possible to set the size of a CSV Layer marker using a field value through the renderer?

Using this EXAMPLE in their samples, they add a CSV layer to the map and set the marker as a simple circle and set the color to Orange using the renderer, but they didn't interact with any field values other than adding them to the Popup. I know it is possible to set Symbol Size on feature Layers, as this example here using a Feature Layer sets the size, rotation and color of an arrow using some attribute values.

I am trying to use their CSV Layer Example to set the size of the Marker by the field 'magNst' here is my attempt:

<!DOCTYPE html>

    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
    <meta name="viewport" content="initial-scale=1, maximum-scale=1,user-scalable=no"/>
    <title>CSVLayer sample</title>
    <link rel="stylesheet" href="https://js.arcgis.com/3.29/esri/css/esri.css">

Trying to use the renderer to access CSV Field Values  doesn't appear to have any affect, am I missing something or is this not possible with the CSVLayer method?

I believe I can try to read/parse the CSV into a FeatureCollection, and the Renderer will work as it does in the example where they rotate the feature layer symbols and change its size. However the ability to add the CSVLayer is comparatively lightweight so I'd like to work with it as much as possible.


texts: Lizmap integration in existing page build with Smarty template

texts:
We start working with Lizmap an are impressed by its capabilities. Since we want to integrate Lizmap in our existing infrastructure, we are looking for a way to add the Lizmap template to our existing Smarty-templates.
Until now we started with a simple iframe structure, which works, but is not with what we would like to end up.
Is there is a way to integrate Lizmap directly in the code of our template?


texts: Sites sector plots on Google Earth using Python appearing elliptical instead of circular

texts:
I am using below function to create sectors of sites for a telecom operator. I am using lines to show different technologies deployed to keep visuals simple. Problem I am having is that sector appears slightly skewed towards elliptical side. Showing as in pictures.

I am feeding out put of this function into SimpleKML which creates KML for my purpose. If any one can guide me to make them proper circular rather than elliptical it will be great.

I have seen a similar question about some one had similar problem while creating a buffer circle around a point in GE. In one of the answers some one resolved it by using PROJ library. From that solution I can't understand how to solve my problem.

def get_points(lon, lat,ht, sectRadius_km, azimuth, shapeType,offsetAzimuth_deg):
earthCircumference_km=40007.863 #40075
sectorRadius_deg = 360*(sectRadius_km/earthCircumference_km)
if shapeType == "Point":
    sectorPointLon_deg = lon + sectorRadius_deg * math.sin((azimuth+offsetAzimuth_deg)*math.pi/180)
    sectorPointLat_deg = lat + sectorRadius_deg * math.cos((azimuth+offsetAzimuth_deg)*math.pi/180)
    #print([sectorPointLon_deg , sectorPointLat_deg])
    return [(sectorPointLon_deg , sectorPointLat_deg, ht)]
elif shapeType == "Line":
    sectorPointLon_deg = lon + sectorRadius_deg * math.sin((azimuth+offsetAzimuth_deg)*math.pi/180)
    sectorPointLat_deg = lat + sectorRadius_deg * math.cos((azimuth+offsetAzimuth_deg)*math.pi/180)
    #print([sectorPointLon_deg , sectorPointLat_deg])
    return [(lon, lat,ht),(sectorPointLon_deg , sectorPointLat_deg,ht)]


texts: Getting latlong in popup on mouse click over linear feature using Leaflet?

texts:
I´ve created a Geojson layer of linear features using the option "onEachFeature" to binding a popup wich shows some data for each feature.
This is the code that is working right for this purpose:

function onEachFeature (feature, layer) {
    layer.id = feature.properties.num,
    layer.bindPopup(
        "
"+
            "Orden de trabajo (" + feature.properties.estado +")
"+
            "
"+
            "Identificador: " + feature.properties.num +"
" +
            "Fecha ini.: " + feature.properties.f_ini +"
" +
            "Fecha est. fin: " + feature.properties.f_estimada+"
"+
        "
"

    )
};

var l_encurso = L.geoJson(<%- ordenesl %>, {
    filter: function(feature, layer) {
        return feature.properties.estado == "En curso";
    },
    style: function(feature) {
        return {
            weight: 5,
            color: 'green'
        };
    },
    onEachFeature: onEachFeature
});

Now, i would like to get the coordinates (lat long) of the mouse position when i click over the linear feature and show them also in the poup.

Can you give me any idea, please?


texts: Merging sections of lines with same geometry in QGIS?

texts:
I have a series of bus routes with the frequency/hour attached in QGIS.

The idea is to combine the frequencies on routes with more than one bus to have an overall bus frequency for the sections of road the buses share.

Is there a method to do this short of manually selecting routes?


texts:
Has anyone ever attempted to create a raster2pgsql gui (like the shp2pgsql gui).

As I am sometimes lazy to use command line, this could be a very nice to have.


texts: Creating thematic map where individual points are basically grid in MapInfo?

texts:
I want to plot a heat map which shows the different usage in different areas. But the thing is, in my original dataset, the data got truncated, which means that each point represents a small area rather than a certain point of latitude and longitude. For example,

(14.21, 55.38)-> represents a square of place for (14.21XXXXXX,55.38XXXXXX)

I have searched the tutorial file and tried, but I can only get my points with different shapes and darkness, rather than making them continuously showed in a map. I am new to GIS and MapInfo.

Can you give me some ideas of how to transfer my points to grids, or how to directly use the points to create heatmap which is not displayed discretely?

A snapshot of my picture:


texts: Attachment match table attaching too many files

texts:
I'm having a problem correctly generating a match file for a folder of pdfs. The pdfs all have numeric file names - 10.pdf, 58.pdf, 586.pdf, 587.pdf, etc. The feature class where the files are getting attached has a field, TxtSt, that contains the name of the pdfs without the .pdf extension - 10, 58, 586, 587, etc. Not every pdf will necessarily match with a value in TxtSt.

When I run the Generate Match Table tool, it creates the match table. The problem is that if the value of TxtSt is 58, the tool is also matching 586, 587, and so on. It's as if the match process uses a wildcard - 58*, for example - to match anything beginning with that value.

I don't know if this is an error on my part or in the tool, but this is puzzling me.


texts: Want to use measure distance and add points feature together for my project

texts: I am beginner in Qgis. I need advise regarding using add point feature and ruler together. I want to draw 10 points each at exact 20m distance on my polygon for animal control program. By using measure line feature (Ruler) I can draw points exact at 20m distance but while i start using add point feature the ruler points are disappearing. how can i fix that? I have QGIS 3.10.6 version. Please help me. Thanks in advance.


texts: Randomly splitting Feature Collection in Google Earth Engine?

texts:
I am trying to do a supervised classification in Google Earth Engine and created around 50 polygons for each of my 5 categories and put them into Feature Collections and merged them into one single feature collection with the property "class".

Now I want to split this feature collection into two random samples - the first sample I will use for training my classifier, and the second for validating the classification afterwards.

I have been trying different strategies. I e.g. tried reducing it to an image (images are able to be split into random samples according to the GEE tutorials), but I could not manage to keep the property "class" (the only relevant property my feature collection has). I also tried inserting a random column in my feature collection and filtering by it, which did NOT work:

//merging feature collections
var fc = water.merge(croplandVEG).merge(croplandBARE).merge(urban).merge(naturalVEG);

var newfc = fc.randomColumn();

print(newfc);

//split the image of the training points into Training & Accuracy samples!

var filterrandom0 = ee.Filter.equals('random',0);
var trainingsample = newfc.filter(filterrandom0);

var filterrandom1 = ee.Filter.equals('random',1);
var accuracysample = newfc.filter(filterrandom1);

// Get the values for all pixels in each polygon in the training.
var training = clipped.sampleRegions({
  // Get the sample from the polygons FeatureCollection.
  collection: trainingsample,
  // Keep this list of properties from the polygons.
  properties: ['class'],
  // Set the scale to get Landsat pixels in the polygons.
  scale: 100
});

// Create a random forest classifier with custom parameters.
var classifier = ee.Classifier.randomForest().train({
  features: training,
  classProperty: 'class',
  inputProperties: bands
});

// Train the classifier.
var trained = classifier.train(training, 'class', bands);

// Classify the image.
var classified = clipped.classify(trained);

// Classify the validation data.
var validated = accuracysample.classify(classifier);

For this I get the error code:

  ConfusionMatrix (Error)

  No data was found in classifier training
  input.

Why is this possible for images but not for feature collections on Google Earth Engine? Here is how I tried reducing it to an image, which also did not work so far:

var newimage = fc
  .reduceToImage(['class'],
  ee.Reducer.firstNonNull()).select(['first'], ['class']);

  // Sample the input imagery to get a FeatureCollection of training data.
var trainingsample = newimage.stratifiedSample({
  scale: 100,
  numPoints: 5000,
  seed: 0
});

var accuracysample = newimage.stratifiedSample({
  scale: 100,
  numPoints: 5000,
  seed: 1
});

To which I get the error code:

  Unable to transform geometry into requested projection.