texts: Highlighting polygons of Mapbox tileset on hover?

texts:
I am trying to highlight polygons from a Mapbox tileset on a map on hover. Based on Mapbox's tutorial I am adding a feature state hover to the data that is set through mousemove and mouseleave events. The data's paint attributes depend on the hover state.

That works fine with a GEOJSON source as in their example. My source is a tileset though, and that does not seem to work properly: polygons that are entirely within one tile are rendered just as expected. Polygons that span several tiles are rendered incorrectly: they change their appearance only until the edge of the current tile. Zooming in and out -- which changes the size of the tiles -- verifies that.

Is it possible to add such a hover effect on data from a tileset?


texts: how can I get the rest attributes like ATTV, ATTL etc. using S57 Reader?

texts:
I'm using GDAL to read S57 files to display them with S52 styles, but I was unable to get the attributions like "ATTV", "ATTL"? Do you know how to get it?

I have put the two supporting files, s57objectclasses.csv, and s57attributes.csv to the correct place following the post http://osgeo-org.1560.x6.nabble.com/OGR-to-access-s57-file-missing-attributes-issue-td3742926.html#a3742929 and get the following fields:

*****************************
OGRFeature(WRECKS):2163
   RCID (Integer) = 2164
   PRIM (Integer) = 1
   GRUP (Integer) = 2
   OBJL (Integer) = 159
   RVER (Integer) = 1
   AGEN (Integer) = 540
   FIDN (Integer) = 2135161787
   FIDS (Integer) = 688
   LNAM (String) = 021C7F43FBBB02B0
   LNAM_REFS (StringList) = (null)
   FFPT_RIND (IntegerList) = (null)

--------------------------------------
   CATWRK (Integer) = 1
   CONRAD (Integer) = (null)
   CONVIS (Integer) = (null)
   EXPSOU (Integer) = (null)
   HEIGHT (Real) = (null)
   NOBJNM (String) = (null)
   OBJNAM (String) = (null)
   QUASOU (String) = 2
   SOUACC (Real) = (null)
   STATUS (String) = (null)
   TECSOU (String) = (null)
   VALSOU (Real) = (null)
   VERACC (Real) = (null)
   VERDAT (Integer) = (null)
   VERLEN (Real) = (null)
   WATLEV (Integer) = 3
   INFORM (String) = (null)
   NINFOM (String) = (null)
   NTXTDS (String) = (null)
   SCAMAX (Integer) = (null)
   SCAMIN (Integer) = (null)
   TXTDSC (String) = (null)
   RECDAT (String) = (null)
   RECIND (String) = (null)
   SORDAT (String) = (null)
   SORIND (String) = (null)
************************************

the only missed is "ATTV", "ATTL", do you know how to get it? am I missing something here?


texts: Data Export to Trimble GPS-> "File Not Usable"?

texts:
So, I've been tasked with updating the features and objects of a large amount of data collected a few years ago in a Trimble GPS. The raw data is no longer in the GPS, and was not backed up anywhere else. However, I do have the original data files (imported into GPS Pathfinder from the Trimble as .SSF files) . But, when I tried to transfer those back to the Trimble, the files in the GPS device as "Not Usable". I looked at the file folders and the .car files for each data set are missing.

How can I make these files usable so I can edit their features, or is there another way for me to access the raw data besides exporting the .SSF files?


texts:
I'm finding difficulties displaying values of a related table in a popup window.

Esri said that at version 3.10 and above, popup templates can include fields from related tables using the syntax (relationships/relationship_id/field_name). I tried this syntax but it is not working.

I also checked the sample in this page: https://developers.arcgis.com/javascript/jssamples/popup_relatedrecords.html

It uses the same syntax to call a related field but it only works under mediaInfos, which is strange!


texts: layer.extrusion to set max height setting Language as VBScript instead of Arcade

texts:
So first I'll show you what I manually set and am now trying to create purely through Python.

I accomplish this by setting the MAX_HEIGHT as shown here

Now this is the cell in my notebook I am using to assign youthHeight ect to my 3D layers.
aprx = arcpy.mp.ArcGISProject("current")
for m in aprx.listMaps() :
    for layer in m.listLayers() :
            if layer.name == "youthCounties3D" :
                layer.extrusion("MAX_HEIGHT","$feature['Counts.youthHeight']")
            if layer.name == "youngAdultCounties3D" :
                layer.extrusion("MAX_HEIGHT","$feature['Counts.youngAdultHeight']")
            if layer.name == "olderAdultCounties3D" :
                layer.extrusion("MAX_HEIGHT","$feature['Counts.olderAdultHeight']")
            if layer.name == "seniorAdultCounties3D" :
                layer.extrusion("MAX_HEIGHT","$feature['Counts.seniorHeight']")

However, when I do this, it comes out a bit different.

The Language is set to VBScript and not Arcade. If I manually switch back to arcade, the expression changes but if I rewrite it properly it works as intended. So it seems to me that I need to keep the language set to Arcade but I don't know how to do that in my code. Why is this happening? How can I set the MAX_HEIGHT through code without the Language changing?
I could be wrong about what I presume the issue to be, but this is as far as I got with this.
P.S. If it helps, below is a screenshot of what happens when I manually change the language back to Arcade after running my cell.


texts: QGIS 2.12 LAStools lasmerge error: "no outputs specified"

texts:
I am trying to merge 4 .las files with "lasmerge" within "LAStools" in QGIS 2.12 running Windows 7 64. I keep on geting error "no outputs specified" even thoug i did specify the output folder and filename. LAStools is placed in C:\LAStools.

Can anyone please provide a solution to a problem?


texts: arcpy raster to polygon without merges (or assign autoincrement values to raster cells)

texts:
I have a raster dataset (tiff) and I want to get a shapefile (polygons) that preserves all cells.

Unfortunately, the Raster to Polygon tool automatically merges all adjacent cells with the same values. I want to prevent this.

The only way that I could think of, is to change all cell values to an auto-incrementing value in order to avoid matching values and therefore a cell merge.

How could I do this? Or is there a more elegant way altogether?


texts:
I am trying to open and extract data from a 90MB TIFF file using Python. The code I'm using is the following:

from osgeo import gdal, osr, ogr

def get_value_at_point(rasterfile, pos):
    gdal.UseExceptions()
    gdata = gdal.Open(rasterfile)
    gt = gdata.GetGeoTransform()
    data = gdata.ReadAsArray().astype(np.float)
    gdata = None

    x = int((pos[0] - gt[0])/gt[1])
    y = int((pos[1] - gt[3])/gt[5])
    return data[x, y]

However, I'm getting the following error message:

RuntimeError: This is a BigTIFF file.  BigTIFF is not supported by this version of GDAL and libtiff.

Any hints on how to overcome this? I'm using Anaconda and my GDAL version is 2.3.3 and the libtiff version is 4.0.10


texts: In ArcGis10.2 ModelBuilder: Make Output-Name in ArcGis System Tool "Featur to Raster" depending on selected field for pixel value

texts:
This might be easy, but I could not find a solution:

I am using the system tool "Feature to Raster" in ModelBuilder in ArcGis 10.2.

I want to substitute the name of the output raster, "RasterName", with a variable output name, depending on the selected "Field" for pixel value. It should be something like Raster_%Field%", which results in "Raster_Used_agric_land


texts:
I have an Openlayers application where I can input a centre point, width and height in meters from which I can generate a bounding box to display on the map.
The problem I'm having is calculating the box with the correct length sides.
The input data is as follows:

Latitude: 53.3530563113
Longitude: -6.45723718258
Height: 50 meters
Width: 50 meters

The lats and longs are converted to EPSG:900913 (3857) and the box creation maths performed to output the following coordinates, shown in both EPSG:900913 & EPSG:4326.

Centre Point X/Y: -718816.355096198/7048572.675410685

Bottom Left Point:
 x: -718841.355096198, y: 7048547.675410685
 Long: -6.45746176140103, Lat: 53.35292226393439

Top Left Point:
 x: -718841.355096198, y: 7048597.675410685
 Long: -6.45746176140103, Lat: 53.35319035824405

Top Right Point:
 x: -718791.355096198, y: 7048597.675410685
 Long: -6.45701260375897, Lat: 53.35319035824405

Bottom Right Point:
 x: -718791.355096198, y: 7048547.675410685
 Long: -6.45701260375897, Lat: 53.35292226393439

The distance between any of the EPSG:900913 vertex points, or the length of any of the vertex is 50 meters (as set) but the distance between the Lat Long transforms is closer to 30 meters.

I'm aware of Cartesian and projected coordinate etc. systems but should the difference be this big or should I do my calculations some other way instead of transforming to 900913 and back to 4326, i.e. only in Lats and Longs to try and get the 50 meters correct in the 4326 system.

Any links/info or maths to get this right?


texts: Making colours proportional to value field, continuously and not discretely as with graduated colours classification in ArcGIS Desktop?

texts:
Under symbology in ArcGIS, there's the option to graduate colours given a value field, a colour ramp, and some classification.

How do you make the colour graduation continuous (rather than discrete as with the classification)? Either continuous along an existing colour ramp or by being able to specify colour values directly through functions for each RGB/HSB channel?

It seems like pretty basic functionality to be able to provide more than the maximum of 32 levels of classification, so I hope I'm just missing something. Maybe I'm not using the right terminology because I haven't found any results by searching yet.


texts: How do I use geometry object in processing.runalg()

texts:
I am trying to clip a polygon (large road network) with every row in another polygon (district polygon) using PYQGIS.

layer_adm_level = QgsVectorLayer(base_path + adm_level + ".shp",'adm_level','ogr')
if not layer_adm_level.isValid():
    print "Layer failed to load!"
for adm_level_feature in layer_adm_level.getFeatures():
    ID = adm_level_feature["ID"]
    try:
        processing.runalg("qgis:clip",adm_level_feature.geometry(),road_network,base_path + "clip_network_" + str(ID).replace(".",""))
    except Exception, e:
        print e

The problem is when I supply the adm_level_feature.geometry() parameter to the clipping function. I get an error Error: Wrong parameter value: <qgis.core.QgsGeometry object at 0x7f98a6311938>.
This means it expects a shapefile, any work around where I can give it an object instead?


texts: Openlayers 3 - Create symbol like openlayers 2 -> OpenLayers.Renderer.symbol.anySymbol

texts:
In the Openlayers v2 I could add symbols like:

OpenLayers.Renderer.symbol.pointToIcon = [100, 70, 20, 50, 42, 70, 20, 90];

How can I do it with Openlayers 3 ?

In the openlayers 3 I can use this feature:

var starSymbol = new ol.style.RegularShape({
    points: 4,
    opacity: .5,
    radius: 10,
    radius2: 10 * .5,
    angle: 130,
    fill: new ol.style.Fill({
        color: "blue"
    }),
    stroke: new ol.style.Stroke({
        color: "red",
        width: 1
    })
});

I can make the other figures,but I can't make the following figure.

In addition, I don't use 'Icon' like:

new ol.style.Style({
    image: new ol.style.Icon(/** @type {olx.style.IconOptions} */ ({
        anchor: [0.5, 46],
        anchorXUnits: 'fraction',
        anchorYUnits: 'pixels',
        opacity: 0.75,
        src: 'data/icon.png'
    }))
})

Because the color of icon is fixed and I can't use my needed colors.


texts: Create image (TIFF/PNG) from geodetic (lat/long) data

texts:
Using Java, I need to render a bunch of geodetic data to an image, ideally TIFF or PNG. Specifically this is weather data, with each sample having a latitude, longitude, and amplitude value.

I'd like to avoid reinventing the wheeling by converting polar coordinates (lat, long) to Cartesian coordinates (x, y in BufferedImage or Graphics2D). Furthermore the amplitude value needs to rendered as a color gradient. I've looks at GeoTools and NASA's World Wind, but it's unclear and not obvious if it has this capability.

Is there a Java library that can do this?


texts: How to show a carto map in an ionic 2 app?

texts:
Has anyone shown a carto map in an ionic 2 app? When I used ionic 1 there was no problem, because it is angular js and I could add the library in the index.html:

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

and then use it in a controller.js like that:

var mapa;
      mapa = new L.Map('mapa', {
        center: [39,-6],
        zoom: 9
      })

But the problem is that now I'm using ionic 2, so it's angular 2, typescript. After adding the libraries I try to do the same code above in the .ts file, but it can't find the name L.

Has anyone idea how to show my map? Or how to use a js library in typescript?


texts:
I will soon have to plan where and how our raster data will be stored as we will be upgrading from ArcGIS 10.0 to ArcGIS 10.2. The following questions first need to be addressed before I can make any demands on my IT Department. The rasters which will be stored are mostly Aerial photos and Topographic Maps. Every year we get new datasets which will be added to the collection. Older versions of the maps will be held in the same database. I work in public administration where around 400 people use a browser based mapping client to call the data. Nobody needs write-Permissions for raster data apart from me.

Where? Using ArcGIS Server I presently have the option of either holding the raster in MSSQL-Server using SDE or holding them in a FGDB.
How should they be served? Direct connect or 3-Tiered as an sde service?
Other recommendations? ECW Format or TIFF? Rastercatalogs? Rastermosiacs?  There seem to be many options available.

Could anyone offer some opinions?


texts: How to export vector file in GrassGIS to QGIS

texts: I am using GrassGIS to analyze watershed. I created a stream vector extracted from raster elevation data, how to extract the stream vector file so I can open it using QGIS?


texts: How can I save a modified shapefile in ArcGIS Pro as new shapefile?

texts:
I downloaded a zipped shapefile and modified it in ArcGIS Pro.

How can I save changes as a new shapefile? I see only the option to save as Layer File.

If it is not possible, how I can edit an existing shapefile and have it in the format of shapefile?


texts: Errors when using QGIS Topology Checker with must not overlap rule

texts:
I use the Topology Checker Plugin in QGIS 3 with the rule "must not overlap" for a single PostGIS multipolygon layer.
Can somebody explain the meaning of errors (colored in red) when the Plugin colors a whole polygon as a red area?


texts:
I have some data that has been mapped in ArcGIS Desktop.

I would like to develop something like if I add/change some parameters a new map can be derived so that it can be used for predicting as well.

All this need be done online.

I need this because mine is a health based project and I want it to be made available for public use.

At least I would like to know how it can be published online.


texts: Why do the range of values change when I export geostatistical layer to raster?

texts: I did an interpolation of soil pH using 232 data points in geostatistical analyst with IDW, RBF and ordinary kriging. For all 3 methods, the range of pH values change (become smaller) after I export the geostatistical layers to raster. I understand that has something to do with the way values are calculated per cell in raster, but is there a way I could keep the pH range in the output raster?


texts: Tables and view not visible in QGIS 3.10.2

texts:
Using QGIS 3.10.2 I created a geopackage DB with 2 geographic layers and some tables, later on I edited the same geopackage with the DB Browser for SQLite 3.12 and created a View, add another table and changed a name of one of the first tables ('Material' was changed to 'Poste_Material'). Later on, In QGIS again,  I re-connected to the geopackage and found that all what I had added/changed with the DB Browser is not visible, not the view not the tables. (I reconnected several times)
I have been looking in the internet and have not found at least the reason why the new table is not visible, neither if there is something wrong with the view.
The SQL I used to create the view is this (I borrowed the syntax from spatialite views):
CREATE VIEW 'CABLES' AS
SELECT C.Id_Cable, P.Nombre_Propietario, L.Id_Linea,  L.geometry as geometry
FROM LINEA AS L
LEFT JOIN CABLE AS C ON L.Id_Linea = C.Id_Linea
LEFT JOIN Propietario as P ON P.Id_Propietario = C.Id_Propietario
WHERE geometry IS NOT NULL
ORDER BY L.Id_Linea, C.Id_Cable


texts:
I get an error message:
"Input layer A contains invalid geometries (Feature 33). Unable to complete intersection algorithm" with the intersect tool.

Now I want to have a look at Feature 33. How can I select it? Apparently the ID doesn't correspond to the feature number.


texts: How to determine if feature A and B are located in the same room?

texts:
I have created a map using MapWindows like the following:

A little explanation here, the is a map of an elderly home. The yellow boxes are room and the white rectangles inside the yellow boxes are the beds.

I have a few things that need to be done using QGIS (preferably integrated with R so that I can automate the process):

I need to check whether the beds are "appropriately" placed in the room. For example, in the "lower-left quadrant" room, the "upper-right quadrant" bed is place outside the boundary of the room, this is not appropriate and I need to fix it.
I need to count how many beds are there in one room. Can I do it in GIS/R?
I need to know if a particular bed is in the same room as another particular bed.

I have about 40 maps and a total of >2000 beds, I can't do the 3 things above manually. I wonder if some mechanical means can help.

I am using MapWindows for creating the maps, QGIS for visualization, and R for data processing (checking the dbf files for the naming of the beds, creation of csv files for giving "information" for each bed). I wonder if the 3 tools will be sufficient for the 3 tasks above.

Thanks in advance.


texts: Change from ArcGIS extent box to lat long pair

texts:
I have a copy of some ArcGIS log files with bounding box extent information. This example point should map to Montana.

Extent:-12711495.260857,6062474.786336,-12711267.136963,6062864.358207

I know that the format is xMin, yMin, xMax, yMax for the points. I've taken the points to make a center point.

xCenter = (xMax + xMin)/2.0
yCenter = (yMax + yMin)/2.0

It seems like these points are in meters and I'm having difficulty finding a formula to translate them into lat long point pairs. Is there a direct formula I can use in my scirpt for processing? I don't have a copy of arcpy to use which seems to be the most suggested way to solve the problem.


texts: Getting total feature count for multiple layers using PyQGIS?

texts:
I want the total feature count of all the features in all my layers. Doing it manually in QGIS 2.18.16 is too slow.

I have been trying PyQGIS, but am not good with PyQGIS.

Any tips?


texts: Calculating new values for column based on another column with Field Calculator in ArcMap?

texts:
I have an ArcMap table with about 10000 records. I want to calculate values for an empty column (Result) based on the values in another column (FVal). You can see a small example of the table below. For each ID, there are four records.  I need to calculate the following ratio for each ID: (Row 1 + Row 4) / (Row 2 + Row 4). For example, for ID = 0, "Result" is (34.6 + 24.5)/(18.4 + 44.7) = 0.9366.

How can I do this with Field Calculator by using Python or VB Script in ArcMap?

ID - Order - FVal   -  Result

0   -  1     -   34.6

0   -  2     -   18.4

0   -  3     -   24.5

0   -  4    -    44.7

1   -  1     -   1041.6

1   -  2     -   1254.3

1   -  3    -  1059.2

1   -  4    -  1395.6


texts: Integrate GeoJSON file with openlayers.protocol.http

texts:
I've got a problem with (the title says it) integrating a geojson file in my openlayers map. I have to use the openlayers.protocol.http way. I examined lots of examples on the web (and on stackexchange) but I cannot figure out what the problem with my exact code is.
it shows the layer in the layerswitcher, but it doesn't show the features.

I hope someone can find the error. Here's my code (relating geojson):

    WebMapping - Übung 12

<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />

<script type="text/javascript" src="http://openlayers.org/api/OpenLayers.js"></script>

<script type="text/javascript">

    var map;

    function initMap() {

        map = new OpenLayers.Map("openlayers");

        var mapnik = new OpenLayers.Layer.OSM();

        map.addLayer(mapnik);

        var mapCenter = new OpenLayers.LonLat(8.692, 49.412);

        mapCenter.transform(
            new OpenLayers.Projection("EPSG:4326"), // transform from WGS 1984
            new OpenLayers.Projection("EPSG:900913") // to Spherical Mercator Projection
        )

         ZoomLevel 15 setzen
        map.setCenter(mapCenter, 15);

        map.addControl(new OpenLayers.Control.ScaleLine());
        map.addControl(new OpenLayers.Control.OverviewMap());
        map.addControl(new OpenLayers.Control.MousePosition());
        map.addControl(new OpenLayers.Control.LayerSwitcher());
    //  map.addControl(new OpenLayers.Control.SelectFeature(cities,{autoActivate:true, hover:true, onSelect:addPopup, onUnselect:removePopup}));

        var geojsonURL = "wikiCities.geojson";

            cities = new OpenLayers.Layer.Vector("Städte", {

                    styleMap: new OpenLayers.StyleMap({
                        "default": new OpenLayers.Style({
                        pointRadius: 5,
                        fillColor: "#ffcc66",
                        fillOpacity: 0.8,
                        strokeColor: "#cc6633",
                        strokeWidth: 2,
                        strokeOpacity: 0.8 } ),
                        "select": { fillColor: "#8aeeef",
                        strokeColor: "#32a8a9",
                        labelYOffset:13,
                        label:"${name}"} //Text entspricht feature.attributes.name
                    }),

                    strategies: [new OpenLayers.Strategy.Fixed()],

                    protocol: new OpenLayers.Protocol.HTTP({
                        url: geojsonURL,
                        format: new OpenLayers.Format.GeoJSON({
                                  externalProjection: new OpenLayers.Projection("EPSG:4326"),
                                  internalProjection: new OpenLayers.Projection("EPSG:900913")
                                })
                    })
            });

        map.addLayer(cities);

    }

    window.onload = initMap;
</script>

<style type="text/css">
    body {
        font-family: "Helvetica Neue", Helvetica, Arial, sans-serif;
    }
    #header{
        height:3em;
    }

    #wiki {
        position: absolute;
        right: 10px;
        height:91%;
        width: 49%;
        top:3em;

    }

    #openlayers {
        position: absolute;
        width: 49%;
        top: 3em;
        left: 10;
        bottom: 2em;
    }

    #info {
        position: absolute;
        height: 1em;
        left: 10;
        right: 0;
        bottom: 0;

        padding: 0.5em;
        text-align: right;
    }
</style>

    Ihr Browser unterstützt keine Inlineframes
    Wikipedia-Integration in OpenLayers

</div>

<div id="info">

</div>


texts: Creating new column based on parts of strings in another column in QGIS

texts:
I have a column with zoning names (with "subzones" like UV036, UMTs2 etc.) and would like to create a new column with the general zone instead (UV instead of UV036, UMT instead of UMTs2... as you see it is not always the same number of letters or the same structure).
I tried to do this in the field calculator:
CASE
    WHEN "c_zonage_b" LIKE 'UH%%%%' THEN 'UH'
    WHEN "c_zonage_b" LIKE 'UE%%%%' THEN 'UE'
    WHEN "c_zonage_b" LIKE 'N%%%%'  THEN 'N'
    ELSE ''
END

(I also specify which column I want this to be in in the field calculator).
% are wildcards to take into account any possible subzone (sometimes the subzone includes numbers or letters and their number vary, but it never is more than four hence the %%%%).
Unfortunately it doesn't work.


texts: Why doesn't the Field Calculator work in ArcGis 10?

texts:
When I try to calculate a field nothing happens.  I have always been able to just select the field, open the field calculator and type = and the new number. This is frustrating because the other GIS person can do it with no problem. I've tried several different tables, numeric and string fields.
I don't know Python and shouldn't need to use it just to change a field of numbers to a different number.  I have uninstalled ArcGis 10 and reinstalled it and SP1.


texts: Running focal statistics with iterator in ArcGIS ModelBuilder?

texts:
I need to extent the Value of a grid (water mirror elevation) within the nodata area until reaching a specific barrier (the limit of the watershed), in ArcGIS 10.2 for Desktop.

I am using this:

Con(IsNull("grid"),FocalStatistics("grid",NbrRectangle(3,3),"MEAN"),"grid")

So far I have applied it successively to the last output, increasing the size of the rectangle slowly. This is run in a ModelBuilder.

I was getting bored and I have tried to run an raster iterator to apply a focal statistics on each Output.

The model runs it only one time and stops. I was expecting to have to force to stop once it reachs the barrier . And the pre-existing value are also interpolated, which is not the purpose here.


texts: Saving vector layer with another coordinate system

texts:
I want to save vector layer with another coordinate system but i still want to have this layer on the same place.

I created layer on the area with not correct coordinate system and after updating to arcgis it is rotated. I want to create new layer (the same polygons) with correct coordinate system on the same place. There is too much work with creating new layer with new coordinates.

Is there posibility to do that what im asking for?


texts: Export a map and a legend as a PNG with OpenLayers 6

texts: I would like to export a map with a legend in the same PNG using OpenLayers 6. I searched in the OpenLayers documentation but I got just  how to export a map  (without a legend). How can I do this?


texts: Identifying coordinate system of historical map for Arctic Ocean?

texts:
I have a historical map (Arctic Ocean) with unknown CRS system.

In order to compare it with my results, I have to merge them somehow.

Could someone suggest the closest possible coordinate reference system in which the map was made or how I could merge my raster with it?


texts: Refresh Coordinate Reference System in QGIS layers

texts:
I have a large project where several files were loaded with its original CRS, not always the same.
To try to speed up my project, I recently converted all .shp files under the QGIS project directory to a single CRS using GDAL's ogr2ogr. Old files were overwritten.
When I check now the CRS in each layer in QGIS, the original CRS is shown.

Is there a way for QGIS to refresh the CRS of the layers?

I'm working with 2.14.14 version.


texts:
ESRI announced that they are releasing the ArcGIS REST API as open technology.  I've got an ISP running Sql Server 2008 with IIS7.

Does anyone know of a good walk through showing how I can write C# using WCF to expose spatial tables compliant with a particular spec (ideally ESRI's)?

Does anyone know where ESRI's spec is?


texts: Iterating clip in ModelBuilder based on different province/county names using ArcGIS Pro?

texts:
I have a featureclass with all counties/provinces from the Netherlands and I would like to do clips for roughly 8/9 counties based on their names in the feature class.

I have a feeling that I have to use one of the iterators for that, but I am not that really familiar with the iterators, so I was wondering which iterator I should use. I tried to use the "Iterate Feature Selection" and then afterwards the clip, however that didn't work. I filled for the null value the name of the counties/provinces, however I didn't get the clipped cities.

I am currently using ArcGIS Pro 2.3.3

(I know it's fairly easy to do manually (since, it's only 8/9 counties(, however I would like to get more proficient with the ModelBuilder)


texts: Open .asc file inside .tgz.tar archive without decompressing with Python and GDAL

texts:
I am trying to open an .asc file inside a .tgz.tar archive with Python and osgeo.gdal.

Here is the link to the official description that this is possible by using /vsitar/

#source file name
src_filename = "/evapo_p_downloads/grids_germany_daily_evapo_p_200901.tgz/grids_germany_daily_evapo_p_200901.tar/grids_germany_daily_evapo_p_20090101.asc"

#driver
driver = gdal.GetDriverByName('AAIGrid')

#open file with gdal
src_ds = gdal.Open("/vsitar/" + src_filename)

but

print(src_ds)

gives me a None. I try several other descriptions for path like /vsitar//path or /vsitar/path with and without "" but nothing seems to be work. What am I doing wrong here?

I Found this solution for opening successfully a tarball

src_filename = "/vsitar/{grids_germany_daily_evapo_p_200901.tgz}/grids_germany_daily_evapo_p_20090101.asc"
src_ds = gdal.Open(src_filename)
print(src_ds)

<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x000001BEEBEDBEA0> >

So it didn't work when I enter the full path like /path/to.tgz/my.tar/file instead I have to remove the .tar path maybe because gdal knows that in a tarball (.tgz.tar) there is only one .tar file allowed.


texts:
Is there are a library open source that supports WFS-T requests and posts ?

I'm looking for an open source library that supports WFS/and WFS/T that works with open scene graph.


texts: Getting OpenLayers Vector features vs clusters

texts:
I am using OpenLayers Vector with a cluster strategy. When I use Vector.features to get the features returned as a list, it returns an array of clusters and not the array of features. This makes the features unable to be addressed individually. Does anyone know how can I get the features list while using the clustering strategy? Any help is appreciated.

(Eg. If there are 5 features in my KML which are clustered as 1, then I get vector.features.length=1. I need the original features list which will have a length of 5)

Below is the implementation of my Vector:

var maps_layer = new OpenLayers.Layer.Vector("KML", {
                     strategies: [new OpenLayers.Strategy.Cluster()],
                     protocol: new OpenLayers.Protocol.HTTP({
                        url: "kml/maps.kml",
                        format: new OpenLayers.Format.KML({
                           extractStyles: true,
                           extractAttributes: true,
                        })
                    })
                 })

thesis.events.register("loadend", this, function(){
      var my_array = maps_layer.features; //Size of this is = no.of Clusters and not features
});

Here is the documentation: http://dev.openlayers.org/docs/files/OpenLayers/Layer/Vector-js.html


texts: About performance of qgis with big layer or a lot of layer

texts:
does anyone has a feedback on the performance of a big layer
by exemple :
I have a layer with 100 features, each feature contains specifics attributes with dynamic data, the data change each second

My question :
is it better to have this kind of layer or is it better to have 100 small layer with just one feature

is it a problem for qgis to manage a lot of layer ?

update :
here is some more information, the feature i use has 8 parameters ( int, double, string ) each is updated every second, 1 layer for 1 feature

I think the maximum I will have is 50 layers
My computer is a dual xeon, 8go of ram, no database, everything is in memory
the application is coded in c++/qt 4.8.5 with qgis api c++ 2.8.2


texts: Creating a dataset from attributes in QGIS Graphical Modeler

texts:
I need to select just a field from a dataset to compose a model in QGIS 3.4 modeler.
How could I extract this field from the main dataset while creating another one?.

I've tried with Drop field(s) but the dataset is to large to enumerate everyfield and make the task efficient.

I also tried Extract by attribute but it returns a blank attribute table with every field from the main dataset but no features.

The output should be similar to a pandas data selection like this:

import pandas as pd
df = pd.read_csv(file_path)
df_select = pd.DataFrame(df['field_name'])


texts:
Is there a python version of qTuio?

qTuio is implemented in C++. However, I am currently trying to use pyqgis due to problems with C++ environment (I cannot even create a raster layer).


texts:
I am using XYZ to Lat Long which I am then converting to UTM as I am unable to find XYZ to UTM implementation directly.So,in attached image there is distortion on boundaries which has particular pattern. I believe this behavior is because of converting Lat-Long to UTM and then using UTM co-od for further processing using gdal.

How do I calculate XYZ or TMS Tiles bounds in UTM directly instead of Lat-Long?


texts: ArcGIS webservice returning only one result when using FindAddressCandidates method

texts:
I am using ArcGIS web service method FindAddressCandidates to fetch matching addresses to the provided input. But the result comes with only the top matching address. All other similar address are ignored and not shown in the result.

Some additional info:

Our primary application is based in Siebel and Siebel right now does not support/accept messages coming in JSON format. Also Siebel does not have an inbuilt converter to change the JSON to XML format. So because of this we had to use the Webservice implementation for ArcGIS instead of making a REST call.

So I was able to get the WSDL file from the url http://resources.arcgis.com/en/help/soap/10.1/index.html#//01vp0000008v000000

Before implementing the whole solution in Siebel, we tried accessing it via SOAP UI. I have attached the request and response message when tested using Web service. You can see there is only one response that has come when accessed using Web Service. Also find below the response which we get with a REST call (response in JSON format), it has multiple address suggestions.

SOAP Request

SOAP Response

REST call RESPONSE

{
 "spatialReference": {
  "wkid": 4283,
  "latestWkid": 4283
 },
 "candidates": [
  {
   "address": "50 LONSDALE STREET, MELBOURNE, VIC, 3000",
   "location": {
    "x": 144.97114435763478,
    "y": -37.809456786877405
   },
   "score": 100,
   "attributes": {

},
   "extent": {
    "xmin": 144.96934571354262,
    "ymin": -37.811255430969553,
    "xmax": 144.97294300172695,
    "ymax": -37.807658142785257
   }
  },
  {
   "address": "LONSDALE STREET, MELBOURNE, VIC, 3000",
   "location": {
    "x": 144.96234167907767,
    "y": -37.812256458234671
   },
   "score": 100,
   "attributes": {

},
   "extent": {
    "xmin": 144.96054303498551,
    "ymin": -37.814055102326819,
    "xmax": 144.96414032316983,
    "ymax": -37.810457814142524
   }
  },
  {
   "address": "LONSDALE LANE, MELBOURNE, VIC, 3000",
   "location": {
    "x": 144.95865963166221,
    "y": -37.812736841323428
   },
   "score": 91.599999999999994,
   "attributes": {

},
   "extent": {
    "xmin": 144.95686098757005,
    "ymin": -37.814535485415576,
    "xmax": 144.96045827575438,
    "ymax": -37.81093819723128
   }
  }
 ]
}


texts: Raster Multiplication - Can't Multiply Sequence by Non-Int of Type "Float"

texts:
I am trying to multiply a list of GRID raster datasets using a different value for each raster. The value is in the raster name. Some examples are [wse_1979_5,wse_1980_0,wse_1980_5], where 1979_5, 1980_0, and 1980_5 will become 1979.5, 1980.0, and 1980.5. Every raster currently is classified where VALUE = 1 for all cells, and multiplying by the above values from the names should make the raster datasets equal to those values. So far I have:

arcpy.env.workspace = *path*
from arcpy.sa import *
rasters = arcpy.ListRasters()
for raster in rasters:
    multValue = float(raster[4:-2]+"."+raster[-1])
    outRaster = raster * multValue
    outRaster.save(raster+"_r")
TypeError: can't multiply sequence by non-int of type 'float'

I am getting "TypeError: can't multiply sequence by non-int of type 'float'". I believe I need to somehow iterate the raster since it is apparently a string(?), but I'm not positive and can't quite figure this one out. It is simple to do in raster calculator, but I have a lot of datasets to multiply.

Thanks!


texts: Publishing Geoprocessing Service that returns error (or without default parameters)?

texts:
I want to publish a GP Service to ArcGIS Server without default input parameter. However, I also want that GP Service to return an error when no input parameter is provided.

How could I do so? ArcGIS wont let met publish a GP Service that returns an error, and if I run it with a valid input value, it will be defined as the default value and I wont be able to remove it.


texts: extend extent of geoJSON layer - Leaflet

texts:
I am zooming to the extent of a set of geoJSONs collected in a FeatureGroup with this call:

 map.fitBounds(myFeatureGroup.getBounds());

...I would like to extend the extent a bit; as it is the extent is right on the edge of the geoJSON.

Is there a way to do this? I am new to Leaflet (I am more familiar with the ArcGIS JS API, where i would use:

map.setExtent(fullExtent.expand(4))


texts: GEE - Confused about error matrix parameters

texts:
I tried to calculating the error matrix for my classifications map.
I am confused about error matrix parameters.

I split my samples into 2 categories, 70% for training and 30% for validation.

Example code:
var validation = validSample.errorMatrix('q1','q2')

Which parameters do I need to use for q1?

The feature collections (30% sample)? or classified image using 30% sample? I mean, do I have to make a new classification image using 30% sample?

q2: classified image using training sample.


texts:
I am trying to use the Global Wind Atlas for a map about alternative energy:
https://globalwindatlas.info
I have downloaded the 14.5GB geotiff of global wind power density. QGIS reads it perfectly and I can apply a color ramp, but I would like to match the existing product exactly. Does anyone know where to obtain a color map or SLD file? Or how to extract the styling information from the online resource?


texts: Matching streets in two line datasets with PostGIS or QGIS?

texts:
Explanation:
I have 2 tables on PostgreSQL, one table named 'city' that have linestrings wich represents streets from a city (using QGIS I downloaded OpenStreetMaps data and imported to PostgreSQL using shp2pgsql).

The second table named 'round' also contain linestring wich represent routes from some vehicles (Using a Google API they were drowed).

Using QGIS I start visualazing how both tables look like:

When I made zoom to compare both tables I realized that differed alot in different positions

My question:
¿Is there anyway to automatically update the lines from table 'round' so they can perfectly match with the lines from table 'city'?


texts: How to autoclamp to ground automatically in Google Earth using a dynamic KML created from ArcGIS Server 10

texts: It loads as an absolute everytime so terrain must be turned off to see the lines.  Any idea how to get it so they default to clamped to ground?  Using newest version of GE.


texts: Export label to dxf from postgresql by ogr2ogr

texts:
I need to export label to dxf from postgresql by ogr2ogr command, label information is text field in PG, but it is always unsuccessful, I always lose the text label information.

  Those commands are here:

    ogr2ogr -f "DXF" "demo.dxf" "PG:host=localhost user=postgres dbname=demo pas
    sword=xxxx port=5432" -sql "SELECT geom, 'LABEL(f:\'Times New Roman\',s:12pt,t:\'' || text || '\')' AS OGR_STYLE FROM demo"

    ogr2ogr -select 'text' -f "DXF" "demo.dxf" "PG:host=localhost user=postgres dbname=demo pas
    sword=xxxx port=5432" "demo"

PG Fields are here:

What is wrong with those commands?


texts: Values on top of each bar of a histogram in QGIS

texts:
Is it possible to show the values on top of each bar of a histogram created with QGIS? (I am using QGIS 2.18)

I know this post [Is it possible to label bar charts? addressed a similar question years ago, but it was unanswered, and maybe the program has evolved since, or someone has find a better solution.


texts: AttributeError: 'unicode' object has no attribute 'replaceWorkspacePath'

texts:
I must of written this code using 15 different structures with no luck. I am simply trying to take an old layer 'Windfarms' stored in 'C:\Users\darle\Desktop\MXD_Test\Test.gdb' and replace it with 'New_Wind_Farm' stored in 'C:\Users\darle\Desktop\New_Shapes\new.gdb'.

my code is as followed:

import arcpy, os
for filename in os.listdir(folderPath):
fullpath = os.path.join(folderPath, filename)
if os.path.isfile(fullpath):
    basename, extension = os.path.splitext(fullpath)
    if extension.lower() == ".mxd":
        mxd = arcpy.mapping.MapDocument(fullpath)
        layers = arcpy.mapping.ListLayers(mxd,'WindFarms')
        for lyr in layers:
            layer_source = lyr.dataSource
            layer_source.replaceWorkspacePath(#New shape path, ....
            "FILEGDB_WORKSPACE","New_Wind_Farm")
            mxd.save()
         del mxd

The error I get is in the title. The documentation states that replaceWorkspacePath is for changing both file path and name. What am I missing?


texts:
Given a polyline L and a point P which belongs to this polyline but that also is not necessarily one of the polyline's vertices, I'd like to measure the length from the start of the polyline L till P.

I got P using an identify task with zero tolerance, so I can assume for right that the point belongs to the polyline.

To reach my goals, I create a fake polyline L' that contains the points P and P', where P'=(Xp + Delta, Yp + Delta). Then I use L' to cut L and get two new polylines L1 and L2. Finally I use one geometry service to get the length I want.

The code snippet to create the fake polyline, cut L and measure the length:

private void IdentifyCompleted(object sender, IdentifyEventArgs e)
{
    if (e.IdentifyResults.Count <= 0)
    {
        return;
    }

    var geometryService = new GeometryService(
            @"http://myarcserver.com/ArcGIS/rest/services/Geometry/GeometryServer");
    MapPoint pLine = new MapPoint();
    // Some workaround here: creating and 'arbitrary' line starting with P
    pLine.X = p.X + 0.1;
    pLine.Y = p.X + 0.3;

    PointCollection pointCollection = new PointCollection();
    pointCollection .Add(p);
    pointCollection .Add(pLine);

    Polyline cutter = new Polyline();
    cutter.Paths.Add(pointCollection);
    cutter.SpatialReference = Mapa.SpatialReference;

    // Transform the polilyne returned from the identify task in a list
    List<Graphic> polylineL = new List<Graphic>();
    foreach (var result in e.IdentifyResults)
    {
        polylineL.Add(result);
    }

    geometryService.CutCompleted += geometryService_CutCompleted;
    geometryService..CutAsync(listaResultado, tesoura);
}

private void geometryService_CutCompleted(object sender, CutEventArgs e)
{
    if (e.CutIndexes.Length == 2)
    {
        // Magically figured out what are the first geometry. It
        // is out of the scope of this question
        Polyline first = getFirstGeometry();
        int theLengthIWant = first.Length();
        // We got it ! ! !
    }
}

I guess is too overload execute one identify task and then one cut operation.

So my question is: is there a simpler to do this? I mean, is there a straight forward way to get the length from the start of one polyline till an arbitrary point on it?

[EDITED]

I found a lot of people using IHitTest to do what I want. Unfortunately, this interface belongs to ArcObjects API which is not available to me right now.


texts: creating vector layer from OSM file using OpenLayers

texts:
Using OpenLayers, I am creating a vector layer over a map.
This vector layer is created from an OSM xml file.

My purpose is to do some experiments in Prolog with this data, for example, querying the shortest path (following the declared "ways" in the OSM file) that will visit certain amounts of POIs, given a starting POI.

My problem is that when I load the nodes and ways in my map basically I see a big mess of features, too close to each other and some of them irrelevant.

My first idea was to filter some of these nodes in the original OSM xml file using XSLT, but after doing this filtering the map did not render the vector layer at all anymore. After some experimentation I realized that the problem was that the "way" elements were referencing non-existing node elements (since I filtered them using XSLT).

I could also have deleted the "ways" elements that reference non existing nodes (using XSLT again), but this does not seem to be a good idea since then I would be loosing valuable information for calculating the shortest path and other graph typical problems.

Then I arrived to the conclusion that I do not want to delete certain nodes, but just trying to make some of them non-visible by default in the map.

Then I use XSLT again to change the "visible" attribute to certain nodes. But the problem is that all the nodes in the file are shown in the map anyway, even if they have the attribute "visible" set to false, as in this example:

<osm version="0.6" generator="CGImap 0.0.2">
  <bounds minlat="50.8319000" minlon="4.3355000" maxlat="50.8599000" maxlon="4.3708000"/>
  <node id="125892" lat="50.8500083" lon="4.3553223" user="eMerzh" uid="15399" visible="false" version="6" changeset="9938190" timestamp="2011-11-24T22:05:32Z"/>
...
</osm>

Below the code I am using to create the vector layer from the OSM file:

  <title>My OpenLayers Map</title>


texts: Converting Gauss-Krueger to Lon/Lat in R

texts:
I have sets of Gauss-Krueger coordinates (zone 3) that I'd like to transform to longitude and latitude in R. Here's an example:

library(sp)
coord_gk <- data.frame(GK_R = c(4555277),
                       GK_H = c(5342348))

coordinates(coord_gk) <- c("GK_R", "GK_H")
proj4string(coord_gk) <- CRS("+proj=tmerc +lat_0=0 +lon_0=15 +k=1 +x_0=3500000 +y_0=0 +ellps=krass +units=m +no_defs")
coord_lonlat <- spTransform(coord_gk, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

> coord_lonlat
SpatialPoints:
         GK_R     GK_H
[1,] 28.97989 47.35438

However, the output is not correct. I expect it to be 12.742421, 48.216707.

Where did I go wrong?


texts: How do I edit a shapefile DBF to repair an extra row?

texts:
Something got screwed up and now the rows of my DBF beyond #1210 are offset by one, meaning that each shape has the attributes that should be assigned to the next shape. If I were in a spreadsheet program I would be able to select rows 1210+, copy, then select 1211+ and paste. But I can't do this in the QGIS attribute table editor, and I can't find an editor that will save shapefile DBFs without damage.

Apparently I have too little mojo to answer my own question this early.

Never mind - I fixed it in Excel. I had previously tried swapping columns in Excel and it didn't work, but I guess rows can be edited.

The only minor issue is that the column titles are now all caps. But all the expressions that call them still work, so I guess that's not a problem.


texts: Make processing.runalg use only selected features

texts:
I have a feature layer and a buffer of that feature layer and am trying to cut out the features from the buffers. This would be trivial, except that the buffers overlay multiple features, so I need to select each feature from the two layers individually, use the difference tool to export them to individual layers, and finally merge them all at the end.

I've written a Python script to automate the selection/difference tool part of this but it seems to be using all features from the layers, and not only the selected ones. I have tested that the code for the selection is working so I believe it is the call of the difference tool that is not working as expected. I understood from the documentation that it should only use the selected features, but perhaps either my understanding or my code is wrong.

I am using QGIS 2.18.3.

import processing

buffer_layer = QgsVectorLayer("buffer.shp", "buffer_layer", "ogr")
patch_layer = QgsVectorLayer("patches.shp", "patch_layer", "ogr")

for feature in buffer_layer.getFeatures():
    buffer_layer.setSelectedFeatures([feature.id()])
    site_name = feature.attributes()[0]
    patch_list = []
    for patch in patch_layer.getFeatures():
        if patch.attributes()[1] == site_name:
            patch_list.append(patch.id())
    patch_layer.setSelectedFeatures(patch_list)
processing.runalg("qgis:difference",buffer_layer,patch_layer,0,"output_%s.shp" % site_name)


texts: How to set attribute value of layer with PHP MapScript?

texts:
How can I set attribute values of a layer object with MapScript?

I'm using MapServer version 5.0.3, the php_mapscript.so lib is loaded in my app context

I have a simple Mapfile:

MAP
    SIZE 500 500
    EXTENT -180 -90 180 90
    UNITS dd
    CONFIG "CPL_DEBUG" "ON"
    CONFIG      "MS_ERRORFILE" "/var/www/Data/Data/mapserver/mapserver.log"
    SYMBOLSET "/var/www/Data/Data/mapserver/symbols.sym"
    FONTSET "/var/www/Data/Data/mapserver/fontset.txt"

#--------------------------------------------------------------------------------
# load lines
#--------------------------------------------------------------------------------

    LAYER
        NAME "loadlines"
        TYPE POLYGON
        DATA "/var/www/Data/Data/mapserver/data/loadlines/loadlines.shp"
        STATUS ON
        OPACITY 100
        DEBUG 5

        CLASS #default coloring
            STYLE
              COLOR 0 0 0
              WIDTH 6
              #PATTERN 40 10 END
            END
            NAME "gust"
            OUTLINECOLOR 254 254 254
            COLOR 100 100 100
            #OPACITY 50
        END

        TEMPLATE "blank.html"
        DUMP TRUE # allow GML export
        PROJECTION
            "init=epsg:4326"
        END
        METADATA
            "wms_title"     "ffg_test"
            "wms_srs"       "EPSG:54004 EPSG:4269 EPSG:4326"
        END
    END
END

Which I want to view by calling this PHP script:

<?php
$map_path="/var/www/Data/Data/mapserver/routeguardTemp/test/";

$map = ms_newMapObj($map_path."loadlines.map");
$layer = $map->getLayerByName('loadlines');
$layer->set('opacity', 10);
$image=$map->draw();
$image_url=$image->saveWebImage();

?>

 <HTML>
  <HEAD>
      <TITLE>Map</TITLE>
  </HEAD>
  <BODY>
      <IMG SRC=<?php echo $image_url; ?> >
  </BODY>
 </HTML>

The script should alter the opacity value of the layer, but there is no effect.

How can I do this?


texts: Selecting features by attribute if in Python list?

texts:
I am trying to complete a select by attribute in Python but based on the query of whether an attribute is present in a list.

Such a query at its simplest should be something like this:

qry = " \"OBJECTID\" in oid_list"
arcpy.SelectLayersByAttribute_management(inft, "NEW_SELECTION", qry)

but that approach returns an invalid expression error.

In the past, I've had to use more complicated sytax for this type of query, such as:

sqlQuery2 = "nid in (" + ','.join(["'"+x+"'" for x in delta_list]) +")"

but an adaptation of this snippet doesn't seem to work for me either, ie.:

 "OBJECTID_1 in (" + ','.join(["'"+str(x)+"'" for x in oid_list]) +")"

What am I missing here?


texts: Ungenerate command in ArcGIS for Desktop?

texts: I am searching for the ungenerate command in argis10 desktop. I get to know that is available only for workstation. But i don't have workstation. Is there any way to use this command or any other way to get the same result as that of ungenerate command.


texts:
I am able to extract geometries from a shape created by the Google Maps DrawingManager class and also put polygons on the map using the the data.loadGeoJson call.

Now I want to do a spatial selection of the polygons that are within the DrawingManager object so I can extract their attributes. Is it possible with the Google Maps API? I know that I can do this by sending the geometries of the shape to a PHP program that can do a query on a spatial database, but I would like to do this completely on the browser.


texts: Error getting GeoJSON data from GeoServer with CQL filter

texts:
I have a problem obtaining GeoJSON data from GeoServer with CQL filter.
When I do the normal query, I get some results, then it's OK ("totalFeatures":2081).
https://mygeoserver/ows?service=WFS&version=1.1.1&request=GetFeature&typeName=test:road&maxFeatures=100&outputFormat=text/javascript&format_options=callback:getJson&SrsName=EPSG:4326&callback=getJson&_=1598536183818

But when I try to use the CQL filter, the query returns a null result ("totalFeatures":0).
My data is in EPSG:4326, comes from an SQL view, and the filter takes as argument the cursor converts to EPSG:3857 to apply a buffer zone in meters.
https://mygeoserver/ows?service=WFS&version=1.1.1&request=GetFeature&typeName=test:road&maxFeatures=100&outputFormat=text/javascript&format_options=callback:getJson&SrsName=EPSG:4326&cql_filter=DWithin(mygeometry,POINT(-66347.340551533%206019052.901573275),500,meters)&callback=getJson&_=1598536183818

Any idea where the problem comes from? Mixing EPSG?

I tried it with a bbox. When I put the parameter directly in the cql filter of the layer on GeoServer, it works; when I call it with leaftlet, it returns an empty result.
cql_filter=bbox(mygeometry,-1.1810302734375002,45.49672163945861,2.6312255859375004,47.12808212013255)

I don't understand where the problem comes from...


texts: How to batch split feature classes using Arcpy?

texts:
I want to bath split analysis tool using Arcpy . I have some feature classes in a feature dataset.I want to Split All of feature classes uisng a shapefile but after split two feature classes, the runtime error raised( Runtime error : Failed to execute. Parameters are not valid. ERROR 000728: Field Province does not exist within table Failed to execute (Split)). My input feature has a province field and the result of the split is fine for a feature class but not for all of them .

How can i solve the problem ?

another qestion is that is there any way to change the output featureclasses? for example province name with underline featureclass name.

import arcpy
arcpy.env.workspace = "D:/split_test.gdb/split"
fcl = arcpy.ListFeatureClasses()
for fc in fcl:
    arcpy.Split_analysis("province",fc,"province","D:/output_split.gdb","")


texts: Using MSID/MID with leaflet to display Google My Maps data

texts:
I'm wondering if there is a way to display MSID/MID with leaflet.

I found this forum below;
Any way I can load KML data from Google maps (again) via JS API?

KMZs work fantastically with the arcgis and google API (kmlLayer) but I can't seem to find any information on using a kmz with leaflet.
I'd prefer not having to Unzip locally, and omnivore does not appear to work with KMZs.


texts: Seeking Great Britain landuse data for free (for academic reasons)?

texts:
For the academic project I am working on, I would need to have the Great Britain landuse data.So far I looked for them in several websites but I did not find anything that was for free and with a good quality (I found some data in osm-x-tractor.org but they are very incomplete).

Is there a good website where I can find these information for free?

Eventually, is there any easy method to extract this large amount of information from OpenStreetMap or any other similar website by myself?


texts: Sorting Relation_Aggregation Expression Results in QGIS3

texts:
I am using a relation aggregation expression to grab a list of inflows from other sewersheds. Some sewersheds take on flows from more sewersheds than others so I approached it through a relation so that I can show all the values from both the Gravity In Flows and the In Flows from other sewersheds. It works but it presents the results in no particular order it seems. Is there something I can add to the expression that will sort the results by the Network Number?

CASE WHEN round((relation_aggregate('Network_Groups','sum',"Network_Consumption     Sheet1_None")),2) IS NOT NULL THEN
'Network '|| "networkGrp" ||'\n   Gravity In Flows '|| round("GFlow",2) ||'     m'||char(179)
||'\n   Liftstation In Flows'||
'\n     '||relation_aggregate('Network_Groups','concatenate', 'Network     '||"In_Network"||' '||"Network_Consumption Sheet1_None"||' m'||char(179),     concatenator:='\n     ')||char(179)||'\n  Total From Liftstations     '||round((relation_aggregate('Network_Groups','sum',"Network_Consumption Sheet1_None")),2)||' m'||char(179)
||'\n  Absolute Total
'||round((relation_aggregate('Network_Groups','sum',"Network_Consumption
Sheet1_None")+ "GFlow"),2) ||' m'||char(179)
ELSE 'Network '|| "networkGrp" ||'\n   Gravity In Flow '|| round("GFlow",2) ||'
m'||char(179) END


texts: How to publish SVG point style to Geoserver from QGIS?

texts:
I used train svg for a point layer (C:\OSGeo4W64\apps\qgis\svg\transport\transport_train_station2.svg) and have used the GeoServer Explorer Plugin to publish the layer on geoserver. But I run into this issue

> Error uploading SVG icon to GeoServer: 404 Client Error: Not Found
> Traceback (most recent call last):   File
> "xxx/.qgis2/python/plugins\geoserverexplorer\gui\explorer.py", line
> 98, in run
>     command(*params)   File "xxx/.qgis2/python/plugins\geoserverexplorer\qgis\catalog.py", line
> 141, in publishStyle
>     self.uploadIcons(icons)   File "xxx/.qgis2/python/plugins\geoserverexplorer\qgis\catalog.py", line
> 157, in uploadIcons
>     raise Exception ("Error uploading SVG icon to GeoServer:\n" + str(e)) Exception: Error uploading SVG icon to GeoServer: 404 Client
> Error: Not Found

I copied SVG file directly to data_dir\styles of geoserver and tried to upload SLD directly (created following burg example) from geoserver interface - still this image is not getting uploaded.


texts: Activating Search Plugin with JS Leaflet

texts:
I am attempting to add the following plugin to search a GeoJSON file I have open in Leaflet Javascript: http://labs.easyblog.it/maps/leaflet-search/examples/geojson-layer.html

I have brought in the search js and css files and called them, I just do not understand the syntax to incorporate it into the existing script.

Here is the script it needs to be able to search (the field is Sheet_Numb):

$.getJSON("Syriashape.json",function(data){
    // add GeoJSON layer and popups to the map once the file is loaded
    L.geoJson(data, {
   style: function (feature) {
       return {};
   },
   onEachFeature: function (feature, layer) {
       layer.bindPopup(feature.properties.Sheet_Numb + "</br>" + "  " + feature.properties.LOC_Catalo);
   }
}).addTo(map);
  });

I know I have to add the following somewhere:

map.addLayer(featuresLayer);
var searchControl = new L.Control.Search({layer: featuresLayer, propertyName: 'name', circleLocation:false});

searchControl.on('search_locationfound', function(e) {

    e.layer.setStyle({fillColor: '#3f0', color: '#0f0'});
    if(e.layer._popup)
        e.layer.openPopup();

}).on('search_collapsed', function(e) {

    featuresLayer.eachLayer(function(layer) {   //restore feature color
        featuresLayer.resetStyle(layer);
    });
});

map.addControl( searchControl );  //inizialize search control


texts: Display together layer in decimal Lat/Lon & layer in DMS

texts:
I have several vector layers with DMS location expressed as something like 5148074,130358 and several vector layers with decimal location expressed as something like 42.23,15.2.

All of them having CRS WGS84 and same is the project CRS.

When I try do display them, they not overlay (although they should): QGIS treats the numbers as is, which means lat 42 and lat 5148074.

How can I fix it?


texts: QGIS Invalid cutline polygon when cutting a raster with mask

texts:
I am working with QGIS 3.10, I should cut some rasters  with a vectorial mask , but whenever I do it, it gives me error: 'Invalid cutline polygon' and the operation provides a raster with pixel values 1.79769e+308, so invalid.
I tried to verify and fix the geometry of the shape file, I tried to do a buffer, and even tried different algorithms in order to do the cut, but the result doesn't change.
Does someone can tell me how to do it?


texts: Processing several KML files in different folders using QGIS

texts:
I currently search for a way to load different KML files into one project. The thing is, that there are several KML files in the folders, which need to be displayed on an OpenStreetMap layer. I have quite less experience using the Modeler or Python scripting and couldn't find any tutorials.
The hierarchy of the folder is like this:
  main_folder
    --> subfolder_1
      --> subsubfolder_1
        --> KML-file ss1_a.kml
        --> KML-file ss1_b.kml
      --> subsubfolder_2
        --> KML-file ss2_a.kml
        --> KML-file ss2_b.kml
      --> KML-file s1_a.kml
      --> KML-file s2_b.kml
    --> subfolder_2
    --> subfolder_3

The folders do hold different amounts of files and each file inside of one folder needs to be visualized in the same file so that I have one result image for each folder. I only want to set the root folder (main_folder) and I want to process the files in each folder recursively.
I was trying to use the different modeler functions, but it is really complex and hard to understand. Could someone give me a hint, how to solve this problem? Is this may be easier to do in Python?


texts: Count valid (unmasked) per-pixel observations in MODIS time series collection

texts:
I want to create a quality image that describes the number of valid (unmasked) observations per pixel time series given a MODIS vegetation image collection.

I want to get a sense for how frequently a given location is masked out because of poor quality or cloud cover to use as a metric to qualify analysis results. For example, if during a given season six observations (images) exist, I want to calculate a map layer that is the count of unmasked observations, so that I can exclude results from pixel time series where the number of observations is less than 3, for instance.

This code block produces a MODIS vegetation index image collection for northern hemisphere summer months and uses the quality band to update each image's mask. How can I count the per-pixel number of unmasked observations given this collection?

// Load all MODIS 13Q1 data from a few months in 2018.
var modisVegCol = ee.ImageCollection("MODIS/006/MOD13Q1")
  .filterDate('2018-05-01', '2018-07-31');

// Define a function to mask image by quality band.
function qualityMask(img) {
  // Make boolean image where high quality pixels are value 1 and all else 0.
  var mask = img.select('SummaryQA').eq(0);
  // Update the image mask and return it.
  return img.updateMask(mask);
}

// Apply quality mask to all images in the collection.
modisVegCol = modisVegCol.map(qualityMask);

// Print the first image.
print(modisVegCol.first());

// Map NDVI mean of the observations.
Map.setCenter(-100.56, 40.97, 4);
Map.addLayer(modisVegCol.select('NDVI').mean(), {min: 0, max: 8000});

Code Editor script


texts: How to use Google Default Pan/Zoom controls in openlayers

texts:
I am using Google maps version 3 with Openlayers,My question is how do you use the default google control for Panning and Zooming and not the default controls Openlayers panning and zooming(The dark blue controls)

I know that it is controlled by the CSS for Openlayers, But what is the default CSS for the Google Pan/Zoom Controls / Where do I find this CSS for the Google controls in the image below.

 var map, layer;

        function init(){

        layer = new OpenLayers.Layer.WMS( "OpenLayers WMS",
                    "http://vmap0.tiles.osgeo.org/wms/vmap0", {layers: 'basic'} );

  // Add Styles
  var myStyles = new OpenLayers.StyleMap({
            "default": new OpenLayers.Style({
                fillColor: "#000",
                fillOpacity: 0.3,
                strokeColor: "#000",
                strokeOpacity: 0.3,
                strokeWidth: 1,
                graphicZIndex: 1
            }),
            "select": new OpenLayers.Style({
                fillColor: "#ccc",
                fillOpacity: 0.6,
                strokeColor: "#ccc",
                strokeWidth: 1,
                label: "Label for Mouseover",
                labelAlign: "cc",
                fontColor: "#333333",
                fontOpacity: 0.9,
                fontFamily: "Arial",
                fontSize: 14,
                graphicZIndex: 2
            })
        });

        vlayer = new OpenLayers.Layer.Vector( "Editable",{
            styleMap: myStyles,
            rendererOptions: {zIndexing: true}
        });

        map = new OpenLayers.Map( 'map', {

        // Add Openlayers Map Controls
        controls:

                 [
                    new OpenLayers.Control.LayerSwitcher(),
                    new OpenLayers.Control.PanZoom(),
                    new OpenLayers.Control.MousePosition(),
                    new OpenLayers.Control.ScaleLine(),
                    new OpenLayers.Control.OverviewMap(),
                    new OpenLayers.Control.EditingToolbar(vlayer)

                 ]
         });

// Make Variables for Google Version 3
var gphy = new OpenLayers.Layer.Google(
    "Google Physical",
    {type: google.maps.MapTypeId.TERRAIN}
);
var gmap = new OpenLayers.Layer.Google(
    "Google Streets", // the default
    {numZoomLevels: 20}
);
var ghyb = new OpenLayers.Layer.Google(
    "Google Hybrid",
    {type: google.maps.MapTypeId.HYBRID, numZoomLevels: 20}
);
var gsat = new OpenLayers.Layer.Google(
    "Google Satellite",
    {type: google.maps.MapTypeId.SATELLITE, numZoomLevels: 22}
);

    // Add layers to map
    map.addLayers([gphy, gmap, ghyb, gsat]);

    map.addLayers([layer, vlayer]);

    // Google.v3 uses EPSG:900913 as projection, so we have to transform our coordinates
    map.setCenter(new OpenLayers.LonLat(10.2, 48.9).transform(
    new OpenLayers.Projection("EPSG:4326"),
    map.getProjectionObject()
), 5);

   var panel = new OpenLayers.Control.Panel();
    panel.addControls([
        new OpenLayers.Control.Button({
            displayClass: "helpButton", trigger: function() {alert('help')}, title: 'Help'
        })
    ]);
     map.addControl(panel);


texts: Getting polygon distance from route using ArcGIS Desktop?

texts:
I need to calculate the location of polygons along the route in ArcGIS (ArcInfo) 9 or 10. What I have already found out is only the calculation of the start and end of the intersection of a polygon with the route. But in my case I need to attribute all the polygons along the route, many of them do not intersect the route at all. The route is one and continuous.

I would need (at best) to get the following values:

Route measure of the polygon starting (aside from the route, by the shortest perpendicular line). Measures to be specified left or right side from the route (as the route measures increase) - probably +/- values
Route measure of the polygon ending (same)
Distance of the nearest polygon edge to the right from a route (zero if intersects)
Distance of the nearest polygon edge to the left from a route (zero if intersects)
Distance of the farthest polygon edge to the right from a route
Distance of the farthest polygon edge to the left from a route

This probably has something to do with VBA, Python or something - if so, please paste the code, since I'm not familiar with any scripting, unfortunately. I bet the problem is common (locating natural hazard sites along a road, for example), but googling and search in ArcGIS help didn't gave me any result. Also, there may be extensions doing that, but my search in ArcScripts didn't yield.


texts: QGIS reports that my .shp file does not have "Protocol" source type

texts:
I tried to add my shapefile, but it reports this error:

Download with "Protocol" source type has failed, do you want to try the "File" source type?

What is wrong with it?


texts: I have downloaded data for india as a dbf file then converted it using osm2po to make it routable but the problem is when i make any routing query it takes a lot of time around 300 seconds and these results are after indexing source, target and id columns. So anyone can please tell how to enhance speed of pgrouting algorithm, or which algorithm should be used, time under 5 seconds would be OK for me, thanks in advance


texts: Using Click Event on Layer A to query Layer B in ArcGIS API for JavaScript?

texts: I have two polygon layers in a webmap service and want to do a dynamic query on click, such that when the user clicks on a lake (layer A) with ID = 100, the corresponding watershed boundary (layer B) with ID = 100 is highlighted. So basically, I want to pass the ID found on click from layer A and use that ID to query layer B by attribute. How do I do this?


texts: "Average" GPS distances between TrackPoints?

texts:
I have a few thousands of GPS tracks (from and for recreational use) and I would like to sort those that were produced by a GPS on the field from those that were drawn (on-line on Google Earth or any other WebGIS or map editor such Garmin Basecamp or others).

Depending on the GPS settings and antenna, trackpoints are usually collected by time or distance to the log file before the final save. On many models this track is usually  compressed (midle trackpoints are erased) to save memory space (I’m pretty sure that some of us still remember that the first affordable models use to have about 20 Mb of memory… enough to take you back or to load a few tracks or routes).

On my data set I have an average of trackpoints distance that goes from 1m to a maximum of 899m (Avg 45.80, Std 50.16).

I have compared this to 396 tracks collected by my Garmin Edge705 on Mountain Biking activities and the average trackpoint distance varied between 10m and 50m (Avg 28.89, Std 7.22). To optimize file size, sport models usually record a new point on every speed or trajectory change making huge variations between trakcpoints distance on the same track, but different models or settings could produce different results.

From my own data, it seems that average distances of up to 50 m could be considered as “collected on the field”, but I believe that GPS APPs for the current smartphones, or other models could push this value a little further…
Within my dataset, tracks with average distance between trackpoints above 200 m are obviously drawn on a medium/small scale map on a WebGIS or other map software. Signs are easily spotted when these tracks are uploaded to a GIS – Tracks don’t usually follow the road, curves are not smooth as if they were collect on the field by a GPS handheld even taking into account GPS errors…

Above 100 m I still “feel” that most are drawn but what could be a reasonable value?


texts: Resetting ObjectID after running Delete_Features on ArcSDE feature table?

texts:
I have a python script that deletes all features on an ArcSDE feature table and then inserts new features to it.

I noted that the ObjectID keeps incrementing from where the last feature was with before deleted.

Just wondering if there is any reset function or tool to make ObjectID start with 1 every time the script runs?


texts: Ogr2ogr driver error when convert shapefile to mif

texts:
I use the following script to convert shapefile to mif:

ogr2ogr -f "MapInfo File" D:\docs\arcgis_data\shp_mif\BOU.mif D:\docs\arcgis_data\shp\BOU.shp

But I got the error:

ERROR 6: Unable to open D:\docs\arcgis_data\shp_mif\BOU.mif

What's the problem?


texts:
I am trying to expose the full path of a file found using os.walk. The walk function is working correctly and is exposing all the files I want however when I try to get the path it looks like it is exposing the path of the current working directory instead of where the file is located.

I tried to get the path by using os.path.join with the original folder name variable but these files are in subfolders

import arcpy, os, time, sys, shutil,codecs
import xml.dom.minidom as DOM
from os.path import basename
from arcpy import env

folderName = "F:\\MXDTest\\"

for root, dirs, files in os.walk(folderName):
    if 'CacheMaps' in dirs:
        dirs.remove('CacheMaps') #don't visit CacheMaps Directory
    for f in files:
        if os.path.splitext(f)[1]==('.mxd'):
            fName = os.path.join(folderName, f)
            print fName

This returns
F:\MXDTest\BaseMapContent.mxd
F:\MXDTest\MunicipalBoundary.mxd

However I want it to return
F:\MXDTest\General\BaseMapContent.mxd
F:\MXDTest\Boundary\MunicipalBoundary.mxd

I want to feed this string into arcpy.mapping.MapDocument()


texts: How to automatically stretch raster in QGIS? (ArcGIS DRA dynamic range adjustment analogue)

texts:
I am looking for ArcGIS's dynamic range adjustment (DRA) analouge in QGIS for raster data sets. And I have found such button:

But I have to click it every time I move raster's extent.

There is any way to automate manual clicking the button?


texts: After how much distance does the x,y coordinate system become noticeably distorted?

texts:
I'm working with maritime sensors that give information of moving objects it detects, in azimuth, range and elevation, and we convert it to x,y coordinates for our application to be able to process it.

We found out that most naval systems consider the earth's center as the point of reference when doing position calculations, but since we needed x,y, we are considering that when our software starts, we note the position of our boat with reference to the earth's center and convert it to x,y coordinates, where the current position of our boat would be 0,0 (the origin).

The boat can move at upto 25 knots, so when we have traveled a good distance, we would convert all x,y coordinates in the software to radius, azimuth and elevation with respect to the earth's center and then convert it to current x,y position of the boat (with the current position being the origin of the x,y grid). Only problem would be, the ambiguity in position of all the moving objects that the sensor was tracking.

Two questions:

Is this the best way to do it or is there a better way? We have to work with x,y coordinates.
How many kilometers would a boat have to move so that the x,y coordinate system would become erroneous due to the earth's curvature?


texts: Build ArcGIS Model to Select by Location (Within a Distance) and Export to KML with an Iteration

texts:
I have a parcel shapefile and point data. I need to select parcels within 0.5 miles of each point and export that data to KML then repeat the set of each point (Iteration).

The KML also needs to be saved as an ID field each point has. I have tried to do this with Model Builder. How would you go about accomplishing this task? I dont know python at the moment but could learn if it is easiest to accomplish this task that way.

This is what I have so far with Model builder:


texts:
I have a tiff image and I would like to georeference it, and so create a Geotiff. To do that, I have to use GDAL. I use in particular, the executable store in the QGIS file (gdalwarp.exe, gdaltranslate.exe...).

For now, I use this code in a DOS command (command line), and this doesn't seem to work.

gdal_translate -of GTiff -gcp 33.9983 669.668 -833589 -692.339 -gcp 124.661 1377.97 -833520 -1379.2 -gcp 6788.34 -10.2993 -826876 6.06375 -gcp 6839.33 748.997 -826813 -755.83 -gcp 3014.52 363.682 -830651 -363.339 -gcp 3892.81 1032.32 -829762 -1021.34 "D:/Users/(...)/Desktop/image_projetee_UTM.tiff" "C:/Users/(...)/AppData/Local/Temp/image_projetee_UTM.tiff"

gdalwarp -r near -order 1 -co COMPRESS=NONE -dstalpha "C:/Users/(...)/AppData/Local/Temp/image_projetee_UTM.tiff" "D:/Users/(...)/Desktop/image_projetee_UTM_georef.tif"

1 - After run this code, the output image is not a GTiff but a Tiff like the input image.

2 - If I try to open the output  image in Google Earth for example, it doesn't seems at all georeferenced

I think I missed something in the way to georeferenced a Tiff file...

Can anyone help me ?


texts: How to make a transparent background image using r.out.png?

texts: Can the white background of the r.out.png of GRASS be changed to transparent?


texts: Can QGIS handle *.adf raster attribute data?

texts:
I have been trying to open Esri grid format (*.adf) files in QGIS. I can display the raster set itself, however, the appropriate attribute data fails to join the spatial dataset.

Can QGIS display the full source raster for ArcGIS?

Or is it a 'proprietary' format that can only be fully displayed in ArcGIS?


texts:
I am trying to compare a set of point features to the centroid of a single feature polygon using ArcPy. The problem is that the polygon generates a set of XY coordinates in meters, whereas the set of points generates a set of XY coordinates in decimal degrees.

Is there anyway to ensure consistency?

Using some variations of SHAPE@XY@DECIMALDEGREES threw errors. When I do the following:

with arcpy.da.SearchCursor("prim_catch_dissolved", ["TR", "SHAPE@XY"]) as cursor:

I get the results in meters, which is not what I'm after.


texts:
Apologies if this question has been answered before.

I'm digitizing a historic coastline and would like to add a single depth figure to it in order to merge it with a current contour model. Is there and easy to add depth to the polyline shapefile I am creating?


texts:
Is it possible to smooth or generalize a polygon in QGIS, perhaps by using a plugin like the generalizer-plugin for lines?

I could convert my polygons to lines, smooth them and convert them back to polygons but this is not working for my polygons because of several reasons.


texts:
Is there any way to represent floors in a building using GIS?

(Also, it's stupid that I have to use existing tags cos none suit my question so I had to choose one at random.)


texts: Add new field in a point layer with an attribute from another layer in QGIS

texts: I'm working with QGIS 3.16 Hannover under MacOS 10.13 environment. I've got a point and a polygon shapefile layers. I want to add a new field to my point layer with the id of the polygon that contains it. Is it possible to do so with the field calculator, instead of creating a new layer?


texts:
I have 30 rasters with 0 and 1 values and I want to combine them into a single raster using the combine tool. However, combine tool has a 20 raster limit, so my here is my current approach and code:

Split the rasters into two batches of 15 rasters and run the combine
Next, run combine on the outputs created in step #1 to create final output.

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

    arcpy.CheckOutExtension("Spatial")

    # Main Directories
    CCDC = "D:/CCDC/"
    CCDC_bi = "D:/LC/CCDC_bi/"

    years = ['1985', '1986', '1987', '1988', '1989', '1990',
        '1991', '1992', '1993', '1994', '1995', '1996',
        '1997', '1998', '1999', '2000', '2001', '2002',
        '2003', '2004', '2005', '2006', '2007', '2008',
        '2009', '2010', '2011', '2012', '2013', '2014',
        '2015']

    for year in years:
        output = os.path.join(CCDC_bi, "year_" + year + "_fix.tif")
        output1 = os.path.join(CCDC_bi, "year_" + year + "_bi.tif")

        if arcpy.Exists(output) and arcpy.Exists(output1):
            print output + " // " + output1 + " >>> Already Exist! "
        else:
            print "Binary Raster Creation Started..."
            outCon4 = Con(IsNull(output), 0, 1)
            arcpy.CopyRaster_management(outCon4, output1,"","","","NONE","NONE","4_BIT","NONE","NONE","TIFF")
            print "Binary Raster Creation Complete"

    comb_list_1 = [] #Placeholder for all combine paths
    comb_list_2 = [] #Placeholder for all combine paths

    bi_years_1 = years[3:-14]
    bi_years_2 = years[-14:-1]

    for i in bi_years_1:
        inRas = os.path.join(CCDC_bi, "year_" + i + "_bi.tif")
        if arcpy.Exists(inRas):
            comb_list_1.append(CCDC_bi + "year_" + i + "_bi.tif")

    for i in bi_years_2:
        inRas = os.path.join(CCDC_bi, "year_" + i + "_bi.tif")
        if arcpy.Exists(inRas):
            comb_list_2.append(CCDC_bi + "year_" + i + "_bi.tif")

    comb1 = ";".join(comb_list_1)
    comb2 = ";".join(comb_list_2)

    start_time = time.time()
    print 'Starting Combine - 1'
    outCombine_1 = Combine(comb1)
    outCombine_1.save("D:/LC/CCDC_bi/combine1.tif")
    print("--- Combine 1 Complete Complete %s seconds ---" % (time.time() - start_time))

    start_time = time.time()
    print 'Starting Combine - 2'
    outCombine_2 = Combine(comb2)
    outCombine_2.save("D:/LC/CCDC_bi/combine2.tif")
    print("--- Combine 2 Complete Complete %s seconds ---" % (time.time() - start_time))

    start_time = time.time()
    print 'Starting Combine Final'
    finalCombine = Combine("D:/LC/CCDC_bi/combine1.tif;D:/LC/CCDC_bi/combine2.tif")
    finalCombine.save("D:/LC/CCDC_bi/combine_final.tif")
    print("--- Combine Final Complete Complete %s seconds ---" % (time.time() - start_time))

I am not sure if this approach will work, I'm currently running this model now to look at the final output. Alternatively, I have tried running combine on these rasters on R but it runs out of memory. Has anyone tried combining large numbers of rasters and if so what method or strategy worked best for them?


texts: Can the QuickOSM QGIS plugin be used with Overpass QL (rather than xml)?

texts:
The QuickOSM plugin for QGIS can support users with the generation of new Overpass queries (for beginners on Overpass like me). It does this using the xml query language. I'd like to also be able to use this plugin with the alternative Overpass QL format. So far I've failed to do this despite having working queries on Overpass-turbo.

Can the plugin support the alternate language. Could I see a couple of queries - written for the plugin specifically - which work so I can base new queries on this?

The main reason for asking is that the answers I've been getting as I learn to use Overpass tend to be in Overpass QL.


texts: Styling Raster with QGIS for publishing in Geoserver

texts:
I have a raster layer with five different landcover types denoting with 1,2,3,4 and 5. in Qgis it looks as follows:

When I tried to upload this layer to geoserver using OpenGeo explorer tools, and visualize it in geoexplorer, the legend changed and showed some unexpectred 'x' there. Why is it so ? How can I retain the same legend style that I made in Qgis  in geoserver? here is the picture from geoexplorer:


texts: trouble restoring dump, $libdir/postgis-2.1 missing

texts:
I have a university project in which I continue a project of last year.
I never used those tool and I'm getting the following error:

ERREUR:  n'a pas pu accéder au fichier « $libdir/postgis-2.1 » : No such file or directory

I download PostgreSQL 9.5 and PostGIS 2.2.1 with stack builder on windows, and try to restore first a "Dump_Dabatase_Structure.sql" and then "Dump_Dabatase_Data.sql" but got the error when trying to restore the structure

EDIT: I deleted postgis with the uninstaller because your command said me i didnt have the postgis extension and used your installer but now im getting:

ERREUR: n'a pas pu charger la bibliothèque « C:/Program Files/PostgreSQL/9.5/lib/postgis-2.1.dll » : unknown error 127
********** Erreur **********
ERREUR: n'a pas pu charger la bibliothèque « C:/Program Files/PostgreSQL/9.5/lib/postgis-2.1.dll » : unknown error 127
État SQL :XX000


texts: Conversion of Latitude/Longitude to different map projection

texts:
I have two m*n matrices coming from a WRF-ARW simulation: one exclusively with longitude values and one with only latitude values. For both matrices I know the originatin projection (lambert conic conformal, with true latitudes and parameter set by me).

I would like to convert them to a regular plate carree projection (the one employed by Google Earth, to be precise): is it possible to do so? If yes, which would be the best instrument and how, exactly?

I'm currently having at hand R (where I've  and may have both QGis and MATLAB (but without any toolbox)

Thanks for any help!


texts: Land Cover Feature Extraction from Satellite Imagery

texts: I'm interested in a low-cost or open source solution for creating land cover GIS layers that utilize both spectral and textural extraction algorithms.  I have used PCI Geomatica, ENVI, and Feature Analyst VLS in the past; however these solutions are a bit beyond my price range, any software recommendations?


texts: Creating layers from selection with while loop using PyQGIS?

texts:
I am trying to write a short snippet of code to break a shapefile up according to a specific attribute.

Here is my code:

from qgis.core import *
import processing
layer = processing.getObject('2012')
#counter is 15 due to limit  R15
counter = 1

# Write the while filter expression and set it
while (counter <= 15):
    selection = layer.getFeatures(QgsFeatureRequest().setFilterExpression('"Group Comp" LIKE  \'R'+str(counter)+'\' '))
    layer.setSelectedFeatures([k.id() for k in selection])

    #Creates layer from selection
    _writer = QgsVectorFileWriter.writeAsVectorFormat(i,r"E:\Scratch\test\NewFile"+str(counter)+".shp","utf-8",None,"ESRI Shapefile", True)
    counter = counter + 1

I get the following error in Qgis when running the above code:

Traceback (most recent call last):
  File "<input>", line 5, in <module>
NameError: name 'i' is not defined

Is my selection of the layer incorrect?


texts: Extracting start and end points from OSM Segment using QGIS

texts: I'm very new to gis. I have osm layer shape file. How can can I extract start and end points of segments and add points id to the attribute table osm segment?


texts: Orientation of linestring with North arrow pointing up

texts:
With PostGIS I am trying to calculate the azimuth of a linestring in order to use it then in a print composer in QGIS, using an atlas. But since the azimuth is calculated by comparing startpoint and endpoint of the linestring, and the startpoint can be further North OR further South than the endpoint, I do not necessarily get a value resulting with the North arrow pointing in the upper 180° of my map.
I would like it to do so, and so my North arrow would never point in any southern direction. I am guessing there is a formula to make absolute the azimuth and so giving me what I need, but I didn't manage to figure that out.
Would anyone have a way to do so, whether it is in QGIS or directly through PostGIS?


texts:
I have lat long coordinates of various merchant ships of the world, say 50000. And I have around 100 instances of coordinates (lat long) of each of these ships at various timestamps.

I have them on my sql db and I wish to plot the tracks of each ship on World Map, not necessarily Google Map.

Please guide how to go about this.


texts: How to remove a pop-up window with PyQGIS

texts:
In my standalone PyQGIS script, at one moment the QGIS CRS window pops-up to select a CRS. I want to add some code to remove the pop-up again when it appears.

I couldn't find anything helpful in the docs yet. Does anyone know how to do that?


texts:
I have Digital Elevation Model (DEM) and elevation points collected from field survey. Now I want to add these elevation points into DEM.

In other words I want to burn these points into DEM. How would I do that in ArcGIS 10.2?


texts: Creating data driven page with two pages per index feature

texts: I want to create a two page 8.5x11 layout for each feature in the index layer. Is there any way to do this in a single DDP? The first page will be created with text from the attribute table and the second page will be a map layout.


texts: What is the difference between a datum and an ellipsoid?

texts:
For example I am converting from ECEF to lla using

def ecef_to_lla(x, y, z):
    lla = pyproj.Proj(proj='latlong', ellps='WGS84', datum='WGS84')
    ecef = pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84')
    lon, lat, alt = pyproj.transform(ecef, lla, x, y, z)
    return lon, lat, alt

I'm not sure how the ellps and datum keywords differ. Documentation points me here: http://proj4.org/parameters.html#parameter-list but it's still unclear how they differ?

I though a datum was a model of earth, usually an elippsoid with different parameters. So what's the difference between a datum and ellps?


texts: Importing symbology for layer from ArcMap to ArcGIS Pro

texts:
I have assigned a symbology to layer in ArcMap and saved the layer as a layer file so that when I open it in ArcGIS Pro, I can visualize using the same symbology but ArcGIS Pro assigns it some random default symbology. I also tried to save it as a layer package file and open it but even that did not work.
How can I move my file between both software while maintaining the symbology constant?


texts: Problem with Cesiumjs sand castle live editor

texts:
I have saved a html file from Cesium SandCastle live editor, but I am not able to run that file on a web browser afterwards, even after loading all the JavaScript files correctly.
How to make it run on the web browser ?
I am using google chrome and webGL is already being installed in Google Chrome.

HTML code-

<!DOCTYPE html>

    <meta charset="utf-8">
    <meta http-equiv="X-UA-Compatible" content="IE=Edge,chrome=1">
    <meta name="viewport" content="width=device-width, initial-scale=1, maximum-scale=1, minimum-scale=1, user-scalable=no">
    <meta name="description" content="Add imagery from a Web Map Service (WMS) server.">
    <meta name="cesium-sandcastle-labels" content="Beginner">
    <title>Cesium Demo</title>

Loading...


texts:
I have run OpenDroneMap(ODM) on a bunch of drone images of a particular area. I have set the Ground Control Points(GCP) using Webodm gcp interface. Now, I have a gcp_list.txt file containing the images with appropriate latitude, longitude values. Now, in the OpenDroneMap Documentaion, it is stated that:

The gcp_list.txt file must be created in the base of your project folder.

Now, I have run ODM keeping the gcp_list.txt file in the images folder. But, I am not clear about where I should put the gcp_list.txt file exactly.

What does this actually mean?

What is the base of my project folder?


texts:
I'm trying to run two Spatialite queries in Python — the first one picks a single point from a table, and the second tries to find all points that are within a given distance of it. So my queries look like:

select geometry from points where point_id = 1

and

select count(*) from points where
PtDistWithin(points.geometry, <geometry from first query>, 100)

I'm stuck on how to inject the geometry that I selected in the first query into the second one. If I drop it in as-is it says cannot concatenate string and buffer objects. I tried selecting it AsText() and there are no errors, but the count is always 0 (there should definitely be matches within 100 meters). Does anyone know a better way of doing this?

Thanks!


texts: QGIS Theme Night Mapping changing text color plugin

texts: Where I have to adjust the color setting in the 'Night Mapping' theme to make the description text from the plugin/extension chooser/updater from black (= unreadable) to white?


texts: Getting RGB values of each pixel in a LANDSAT 8 GeoTIFF with Python 3

texts:
I'm working with a set of LANDSAT 8 Collection 1 Tier 1 GeoTIFFS pulled from Google Earth Engine, and I need to find the RGB values of every pixel in this GeoTIFF.

I've tried using the code from this answer to a similar question:

from osgeo import gdal
gdal.UseExceptions()

ds = gdal.Open(fname)
band = ds.GetRasterBand(1)
arr = band.ReadAsArray()
ct = band.GetColorTable()

# index value to RGB (ignore A)
i2rgb = [ct.GetColorEntry(i)[:3] for i in range(ct.GetCount())]

But I'm getting the following error:

Traceback (most recent call last):
  File "get_color.py", line 12, in <module>
    i2rgb = [ct.GetColorEntry(i)[:3] for i in range(ct.GetCount())]
AttributeError: 'NoneType' object has no attribute 'GetCount'

Ideally, I'd like to have a matrix of (R, G, B) tuples or lists, but I'm open to other formats.

Here's my gdalinfo dump:

Driver: GTiff/GeoTIFF
Files: LC08_014030_20170220.tif
Size is 749, 1570
Coordinate System is:
PROJCS["WGS 84 / UTM zone 18N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-75],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32618"]]
Origin = (604200.000000000000000,4855080.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  604200.000, 4855080.000) ( 73d42'13.61"W, 43d50'29.23"N)
Lower Left  (  604200.000, 4807980.000) ( 73d42'46.30"W, 43d25' 2.88"N)
Upper Right (  626670.000, 4855080.000) ( 73d25'27.68"W, 43d50'16.59"N)
Lower Right (  626670.000, 4807980.000) ( 73d26' 7.41"W, 43d24'50.42"N)
Center      (  615435.000, 4831530.000) ( 73d34' 8.82"W, 43d37'40.10"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
  Description = B1
Band 2 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = B2
Band 3 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = B3
Band 4 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = B4
Band 5 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = B5
Band 6 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = B6
Band 7 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = B7
Band 8 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = B8
Band 9 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = B9
Band 10 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = B10
Band 11 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = B11
Band 12 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = BQA


texts:
I have Polygons with WGS84/EPSG:4326 coordinates
(like 13.57933804767844,52.5126298984133)
and use the Polygon.getArea function to obtain the area of the polygon.
But the result I get is something like 7.695432394893684E-8 .

There were similiar questions (see here, here or here)
but although some of them are marked as answered,
there is not an answer which helped me.

I would like to convert the result to square meter.
How can I do this?


texts: Iterative raster calculation script using NumPy arrays breaks after ~10 rasters, returns: "TypeError: Cannot create raster for numpy array."

texts:
Here is my code:

#Generate Cumulative Thresholds
#5/15/14

import os
import glob
import arcpy
import numpy as np

#Enable overwriting output data
arcpy.env.overwriteOutput = True

#########################################################################################################
def cumulativecalculation(i):

    #Set geoprocessing variables
    inraster = i
    des = arcpy.Describe(inraster)
    sr = des.SpatialReference
    ext = des.Extent
    ll = arcpy.Point(ext.XMin, ext.YMin)

    #Convert GeoTIFF to numpy array
    a = arcpy.RasterToNumPyArray(inraster)

    #Flatten for calculations
    a.flatten()

    #Find unique values, and record their indices to a separate object
    a_unq, a_inv = np.unique(a, return_inverse=True)

    #Count occurences of array indices
    a_cnt = np.bincount(a_inv)

    #Cumulatively sum the unique values multiplied by the number of
    #occurences, arrange sums as initial array
    b = np.cumsum(a_unq * a_cnt)[a_inv]

    #Divide all values by 10 (reverses earlier multiplication done to
    #facilitate accurate translation of ASCII scientific notation
    #values < 1 to array)
    b /= 10

    #Rescale values between 1 and 100
    maxval = np.amax(b)
    b /= maxval
    b *= 100

    #Restore flattened array to shape of initial array
    c = b.reshape(a.shape)

    #Convert the array back to raster format
    outraster = arcpy.NumPyArrayToRaster(c, ll, des.meanCellWidth, des.meanCellHeight)

    #Set output projection to match input
    arcpy.DefineProjection_management(outraster, sr)

    #Set the output filename
    intername1 = os.path.splitext(os.path.basename(i))[0]
    filename = intername1[5:]
    finalname = filename + "_cumulative" + ".tif"

    #Save the raster as a TIFF
    outraster.save("E:\\NSF Project\\Salamander_Data\\New_Cumulative_Rasters\\" + str(finalname))

    #Announce file creation
    print("Created: " + str(finalname))

#########################################################################################################
#Loop through files in directory
geotiffs = glob.glob("E:\\NSF Project\\Salamander_Data\\NoDataToZero\\HadleyGCM\\*tif")

for i in geotiffs:
    cumulativecalculation(i)

This works perfectly for about 10 rasters, then "breaks."  I get the error:

  File "E:/NSF
  Project/Salamander_Data/PythonScripts/Calculate_Cumulative_Raster_Iterative.py",
  line 73, in
      cumulativecalculation(i)

  File "E:/NSF >Project/Salamander_Data/PythonScripts/Calculate_Cumulative_Raster_Iterative.py",
  line 52, in cumulativecalculation
      outraster = arcpy.NumPyArrayToRaster(c, ll, des.meanCellWidth, des.meanCellHeight)

  File "C:\Program Files
  (x86)\ArcGIS\Desktop10.2\arcpy\arcpy__init__.py", line 1860, in
  NumPyArrayToRaster
      return _NumPyArrayToRaster(*args, **kwargs) TypeError: Cannot create raster for numpy array.

I'm not sure why this particular file seems to be causing trouble.  Even if I end up with missing data, I'd at least like to be able to skip over input like this in order to process the rest of the rasters (around 400).

Any ideas on how to do that, or what might be throwing a monkey wrench in arcpy.NumPyArrayToRaster ?