texts: Error using random raster for raster calculations

texts:
I have created a random raster in ArcGIS Pro using the following Python command:

arcpy.CreateRandomRaster_management(Path, "rnraster", "INTEGER", "", 10)

Which works fine.

Performing map algebra operations on this raster is possible when the second parameter is an integer such as:

addition = Plus("rnraster", 6)

However, performing operations that involve any existing raster fail with an enexpected error, such as when I try to add this random raster to my elevation raster/DEM:

addition = Plus("rnraster", "elevdata")
OR
addition = Raster("rnraster") + Raster("elevdata")

I know that the elevdata raster works in map algebra expressions with other rasters, so my question is:

What must I do to the randomly generated raster to enable map algebra operations with existing rasters?

I have already tried matching bit depth, cell size and raster format (16 bit unsigned, 10 and Esri Grid respectively). Does the random raster require spatial referencing first?

Error Log:

    Traceback (most recent call last):
  File "Monte Carlo Simulation.py", line 29, in <module>
    demError = arcpy.Raster("elevdata") + arcpy.Raster(rnRaster)
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\sa\Functions.py", line 4328, in Plus
    in_raster_or_constant2)
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\sa\Utils.py", line 53, in swapper
    result = wrapper(*args, **kwargs)
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\sa\Functions.py", line 4325, in Wrapper
    ["Plus", in_raster_or_constant1, in_raster_or_constant2])
RuntimeError:  ERROR 999998: Unexpected Error.
Failed to execute (MonteCarloSim)

Elevation Data Raster Properties:

Random Raster Properties:

The random raster has no coordinate system or projection set, whereas the elevation raster does, but would this cause the problem?


texts:
I have a SpatialLinesDataFrame and a SpatialPointsDataFrame. I want to calculate the number of points that fall on a line by a specific field in the lines data frame.

Example: I have a river network broken up into 125 fragments with a specific length: I want to count the number of points that fall in each of the 125 fragments. Is there a function to do this like poly.counts?


texts: Checking numeric values in a field and then grouping it out in the where clause

texts:
Building a desktop application using C#.net and arcobjects. I'm trying to retrieve records from a shape file using IQueryFilter. Having trouble writing the where clause in my scenario.

The scenario is as follows.

There are 3 attribute field called "StreetName", "TypeID", "Allowed" in the shape file. The "StreetName" field can have street names such as "Main Avenue", "Robert street", "5th Street" , "ST 6N","CR 7" etc. I need to only retrieve the streets which have a number in it. ( so in this case it will be just "5th Street" , "ST 6N","CR 7") . Once the streets with numbers are retrieved I need to group it according to the "TypeID" field. The "TypeID" field is of type double. Then I need to check for a given "TypeID" whether all the records in the "Allowed" field is the same. ( The "Allowed" field just has 2 different values either Yes or No) If there are records that are not equal in the "Allowed" field for a given "TypeID" I need to get those as the query result.

Example values from the attribute table. ( Record Number column is added for illustration purposes only)

So according to the example data, record 4 ( which has a Yes in the "Allowed" field where all the others of that "TypeID" has No) and record 9 ( which has a No in the "Allowed" field where all the others of that TypeID has Yes) should be the result of the query.

I can select the numbered streets by using a where clause like this.

IQueryFilter spatFilt;

spatFilt.WhereClause = "ST_NAME LIKE '%1%' OR ST_NAME LIKE '%2%' OR ST_NAME LIKE '%3%' OR ST_NAME LIKE '%4%' OR ST_NAME LIKE '%5%' OR ST_NAME LIKE '%6%' OR ST_NAME LIKE '%7%' OR ST_NAME LIKE '%8%' OR ST_NAME LIKE '%9%' ";

After that I'm stuck on how to group according to TypeID and then select the invalid records..

Any help in writing the where clause in the query will be greatly appreciated. (I need only the invalid records. So in this case it should be just 4 and 9)


texts:
Everything is in the title, I would like to save all the tiles from a LAScatalog in one las file. I don't want to store all the tiles in one folder, just one file (because I found it more convenient for the organization of my storage device).
Is there a better way than a dirty loop loading the tiles and merging with rbind.LAS?


texts: How to change visibility when button is clicked again in PyQGIS?

texts:
I am trying to create button that change visibility of grouped layers. When person click button all layers inside group is hidden when is clicked again all layers are visible.

So far i know how to hide layers and set them visible in two separate buttons. But i cannot add to function inside one button.

     #connect to Hide_int button that hide intensity rasters
     self.dockwidget.Hide_int.clicked.connect(self.hideint)

     #connect to Hide_img button that hide img rasters
     self.dockwidget.Hide_img.clicked.connect(self.hideimg)

    def hideint(self):
        self.maps_treeview = QTreeView()
        QgsProject.instance().layerTreeRoot().findGroup('Int').setVisible(0)
        print "Connected"
        pass

    def hideimg(self):
        QgsProject.instance().layerTreeRoot().findGroup('Int').setVisible(2)
        QgsProject.instance().layerTreeRoot().findGroup('Int').updateChildVisibility()
        print "Connected"
        pass

I am using Qgis 2.16. Can you help me?


texts: Severe performance difference using deleteSearchedRows

texts:
The issue I am having involves two versioned tables.  A function deletes 48 rows from Table A and and 48 rows from Table B using deleteSearchedRows.  However, the delete from Table A takes roughly 4 seconds while the delete from table B takes roughly 67 seconds.

The query filter for both deletes is a single number field that has been indexed in the table, objectID indexed as well.  Table A has roughly 45,000 records and Table B 60,000 records.  Both tables also have cascade deletes set to No, so that the only data being deleted would be that table's.

Table structures are very similar.  Table A has OID field + 3 Long integer fields, 2 short integer fields and a date field.  Table B has OID, 3 long integer, 1 Text, and a date field.  Running the delete query from SQL Plus is very quick, so it's probably not a drive efficiency issue.

I broke apart the select and delete sections to see if the search was slow.  I used the Search(queryFilter, false) and it returns instantly, but looping and deleting every IRow.Delete takes slightly over a second per row for Table B.

I am at a loss on what the difference could be to cause such a severe performance difference in two tables.

Using ArcGIS 10.0 with an Oracle database


texts:
I am using r.series command in qgis. The command is performing correctly, but when I import input raster files, it creates a list in a random order. Sometimes the order of input rasters is essential to the outcome. Is there a way to control the list of input rasters?

the order of raster is random


texts:
I'm using this statement in an add in to update the first pane of the ArcMap 10 statusbar:

ArcMap.Application.StatusBar.set_Message(0, message);

When the method is called, no exceptions are raised.  In my OnClick method, I display a message box and when that occurs, I do see the text of the last message.  It seems like the pane needs to be refreshed, but I don't see a method that does that.

I tried:

 ArcMap.Application.RefreshWindow();

but that doesn't help.


texts: What's regular blocking arrangement ? (GDAL)

texts:
I am using GDAL tools to retrieve PostGiS raster. The code(command line) I am using is as follows :

gdal_translate -of AAIGrid "PG:host=localhost dbname='database' user='postgres' password='admin' schema='public' table='rastertable' where='rid=2'" postgis.asc

And it shows error:  rastertable doesn't seem to have regular blocking arrangement. Only tables with regular blocking arrangement can be read from now.

I never heard the blocking arrangement before and I am sure the table in my postgis is a raster table. Any ideas? Thanks.


texts: Preparing an integrated band depth image

texts:
How do I prepare an integrated band depth image?

If I have to prepare an 1 micrometer integrated band depth image (sum of 0.79 and 1.3 band depth micrometer bands),  do I just average all bands from 0.79 and 1.3 micrometers?


texts:
I'm trying to create an interpolation map visualizing the global distribution of genetic markers - basically replicating the attached map. Can it be done in QGIS 3.8 (Zanzibar)?

I installed QGIS 3 but don't see the 'Interpolation plugin' as downloadable (experimental & deprecated plugins enabled). Is it unavailable or found elsewhere?


texts: 3d vector and geometry operations in open-source GIS

texts:
There is some really good open-source software available to work with 2D vectors. For example QGIS has operations like buffer / clip / dissolve / select by location / join attributes by location... But is there something allowing these operations in 3D?

The closest from what I have tried so far seems to be GRASS. It handles 3D points, 3D lines and... Possibly 3D areas as boundaries? There are 3D faces but from GRASS manual it doesn't seem to be supported by v.select etc. And according to http://grasswiki.osgeo.org/wiki/Help_with_3D the volumes doesn't seem to be supported at all at this moment.

For operations with surfaces (3d areas) raster is proven to be reasonably good alternative and are well supported. But for volumes I even can't find a way how to convert closed face mesh / volume to 3D raster + the processing time and data amount.

Specifically is there any open source tool able to do operation "3d points within volume"? (volume defined by faces in DXF format)

Maybe PostGis or some Python libraries able to load volumes?

But furthermore tool able to do volume overlays, buffers, join attributes by location, clip by plane, clip lines / points / areas by volume etc.?


texts: Classification using hyperspectral or Multi-spectral data

texts:
Just want to ask one question regarding the classification using the hyperspectral or Multi-spectral data, If we have hyperspectral or Multi-spectral data, and we only use the data of  R G and B frequencies and we got the good classification result. Does it make sense?
Or we have to use all the possible the frequencies which have the information.


texts: Sort related objects in Identify window of ArcMap 10.2?

texts:
I have a table Points related to a table Events. The event table contains the fields Point-ID, Event type and Year. The two tables are related to each other.

If I choose a number of points with the identify tool and expand the tree in the identify window I can see all the events that are related to a certain point. Is there any way I could make ArcMap sort the related events after year? Right now the list can go 2004, 1973, 1982, 2005 which isn't very pretty.


texts: Calculating radius / buffer on EPSG 3785 in meters vs buffer in WGS 84 in meters

texts:
The goal is working in EPSG 3785 and in particular Google Maps, one wishes to render (on a tile) a circle of a given radius in meters. The rendering is done in pixels.

I would guess that one could take a very simple linear approach to identify the meters per pixel for a given tile zoom:

metersPerPixel = MercWorldWidthMeters / (2 ^ zoom) / tileSize

Now, considering you wish to draw a radius of let's say 100 meters around a point in pixels, then you convert to pixels:

pixelRadius = metersRadius / metersPerPixel

This should in my mind produce a circle of metersRadius on the tile.

However, actually measuring the produced circle's radius with a circle produced in WGS 84 and Great Circle logic, it shows that it is actually about 40% smaller than expected... a big difference.

Does the logic above have a flaw?

Note that I am not interested in finding an alternative algorithm (going to Great Circle calculations and the like)... my question here is "why are the results (so much) different than expected".

All conversions this far have been working fine with the above logic but once I decided to calculate actual distance on EPSG 3785 I seem to be hitting a break wall.

To give an example with figures:

Radius: 0.08 miles ( = 128.74752 meters )
MercWorldWidthMeters: 20037508.342789 * 2
Zoom: 16
TileSize: 256

Above logic suggests a circle of ~107 pixels diameter.
Great Circle suggests a circle of ~175 pixels diameter.


texts: How many GPS-Satellites are needed for time synchronization without positioning

texts:
Is it possible to do time synchronization without positioning with only 1 GPS satellite?
If not, how many satellites do i need for time synchronization without positioning?


texts: Call Vehicle Routing Problem Service via REST - Problem with "time" attributes?

texts:
I am trying to get a Vehicle Routing Problem to work as a geoprocessing service on one of our servers.

If I call the geoprocessing service via REST url and submit depots, orders and routes WITHOUT time information, output stops are returned. As soon as I also submit time information (TimeWindowStart1 etc.), no output stops are returned.

I have tried different time formats: "10:00:00", "10:00 AM", Unix Timestamp (with and without the "), but no success.

Before I publish the service, the VRP works on my ArcGIS Desktop with times.

Our ArcGIS Server configuration: ArcGIS Server 10.1.1 with Security Patch 1, System: Red Hat Enterprise Linux 6.3, 64 bit

Here's one of the REST commands that I tried. Can anyone spot a mistake?

http://arcgisserver01.cloud.alta4.com:6080/arcgis/rest/services/koop_uni_fh_trier/TrierApp/GPServer/SolveVRP/submitJob?
orders={
 "displayFieldName": "POIs",
 "geometryType": "esriGeometryPoint",
 "spatialReference": {
  "wkid": 3857,
  "latestWkid": 3857
 },
 "features":[
 {
    "attributes":{
        "Name":"POI1",
        "TimeWindowStart1": 1367820000,
        "TimeWindowEnd1": 1367848800,
        "MaxViolationTime1":0,
        "ServiceTime":10
    },
    "geometry":{"x":739116,"y":6404262.603}
 },
 {
    "attributes":{
        "Name":"POI2",
        "TimeWindowStart1": 1367820000,
        "TimeWindowEnd1": 1367848800,
        "MaxViolationTime1":0,
        "ServiceTime":10
    },
    "geometry":{"x":739216,"y":6404262.603}
 }
 ]
}
&routes={
 "displayFieldName": "Route",
 "features": [
 {
    "attributes":{
        "Name":"myRoute",
        "StartDepotName":"myStartPoint",
        "EndDepotName":"myEndPoint",
        "EarliestStartTime": 1367820000,
        "LatestStartTime": 1367820000,
        "MaxOrderCount": 30
    }
 }
 ],
 "exceededTransferLimit": false
}
&depots={
 "displayFieldName": "Start- and Endpoint",
 "geometryType": "esriGeometryPoint",
 "spatialReference": {
  "wkid": 3857,
  "latestWkid": 3857
 },
 "features": [
 {
    "attributes":{
        "Name":"myStartPoint"
    },
    "geometry":{"x":739612,"y":6404755}
 },
  {
    "attributes":{
        "Name":"myEndPoint"
    },
    "geometry":{"x":739612,"y":6404755}
 }
 ]
}


texts: Where can I find documentation on the NMEA sentences supported by the Garmin GLO?

texts: I just picked up a Garmin GLO Bluetooth GPS device. Does anyone know where I can find documentation of the NMEA sentences it outputs?


texts: Use of Geoserver raster processes within OpenLayers: WCS/WPS problem

texts:
I've been working with some GeoServer processes using OpenLayers as client. I've been successful in those which use WFS as input (like gs:BufferFeatureCollection). My problem comes when I want to use processes that need as input raster data (tiff). I know that OpenLayers is not a WCS client, but you can send a POST request to the server to get back the result, right?. Particularly I want to execute the gs:RasterZontalStatistics, that needs a tiff and a WFS layer as input, and returns some statistics (average, min, max, etc.) for the tiff inside the geometries of the WFS Layer.

This is the code I've used so far:

/*GEOPROCESS - Calculate the average/min/max values of the MDE raster according to user drawing selection*/
function obtainElevationProcess(){
    var wpsFormat= new OpenLayers.Format.WPSExecute();
    var doc= wpsFormat.write({
        identifier: "gs:RasterZonalStatistics",
        dataInputs:[{
            identifier:'data',
            reference:{
                mimeType:"image/tiff",
                href: "http://geoserver/wcs",
                method: "POST",
                body: {
                    wcs: {
                        identifier:'forestry:mde200_V4',
                        version:"1.1.1",
                        domainSubset: {
                            boundingBox: {
                                projection: 'http://www.opengis.net/gml/srs/epsg.xml#3857',
                                bounds: new OpenLayers.Bounds(-200931.8344047, 4681799.186142525,32948.35238090271, 4948113.428142309)
                            }
                            },
                        output: {format:'image/tiff'},
                        }
                        }
                        }
                        },
            {
                identifier:"zones",
                reference:{
                            mimeType:"text/xml; subtype=wfs-collection/1.0",
                            href: "http://geoserver/wfs",
                            method: 'POST',
                            body:{
                                wfs:{
                                    version:'1.0.0',
                                    outputFormat: "GML2",
                                    featureNS: "http://localhost:8080/thesis/forestry",
                                    featureType:'forestry:Lics',
                                }
                            }
                        }
            }],
            responseForm:{
                    rawDataOutput:{
                        mimeType:"text/xml;subtype=wfs-collection/1.0",
                        identifier:"statistics"
                }
            }
            });

    OpenLayers.Request.POST({
        url:"http://localhost:8080/geoserver/ows/",
        data:doc,
        headers: {
            "Content-Type": "text/xml;charset=utf-8"
        },
        success: function(response){
            alert(response.responseText);
            }
    });

}

Once executed, OpenLayers returns this message (selection):

Process execution ecacba6c-a676-46b5-9010-bd4a2bb35632 failed
Failed to retrieve value for input data
Unrecognised request type {service=WCS, Output=image/tiff,   Identifier=net.opengis.ows11.impl.CodeTypeImpl@1ad2fc (value: forestry:mde200_V4, codeSpace: null), DomainSubset=net.opengis.ows11.impl.BoundingBoxTypeImpl@111d0b5 (lowerCorner: [-200931.8344047, 4681799.1861425], upperCorner: [32948.352380903, 4948113.4281423], crs: http://www.opengis.net/gml/srs/epsg.xml#3857, dimensions: null), version=1.1.1}</ows:ExceptionText></ows:Exception></ows:ExceptionReport></wps:ProcessFailed></wps:Status></wps:ExecuteResponse>

Does anyone an idea about how this can be done (if possible)?

My environment is: GeoServer 2.3.5 / OpenLayers 2.13 / localhost on Apache.


texts: OpenLayers map breaks at certain zoom level when many popups exist in the map

texts:
I'm still learning OpenLayers. I found this weird issue in my OpenLayers map. I create FramedCloud popups and add them to each features.

popup = new OpenLayers.Popup.FramedCloud("popup",
    OpenLayers.LonLat.fromString(parsedFeature.geometry.toShortString()),
    null,
    "
date: " + someDate + "
longitude: " + data.lon + "
latitude: " + data.lat + "
",
    null,
    false
);

But I don't add them until it gets to certain zoom level. Then I add the popups to the map with the following code

olMap.events.register('zoomend', olMap, function () {
    var zoom = olMap.getZoom();
    var popup;
    if (zoom > 11) {
        for (i = 0; i < dm.currentLocationsVector.features.length; i++) {
            popup = dm.currentLocationsVector.features[i].popup;
            if (!popup.map) {
                olMap.addPopup(popup);
            }
        }
    } else if (zoom < 12) {
        for (i = 0; i < vector.features.length; i++) {
            popup = vector.features[i].popup;
            if (popup.map) {
                olMap.removePopup(popup);
            }
        }
    }
});

Actually it works good until zoom 17 then when the map zoom level gets 18 the map breaks.

It looks like the popups and the map get messed in some way.
I don't think I have problem in creating popups.Thought it might some bug in OpenLayers.


texts: How can I add another WHERE filter to a PostGIS query?

texts:
I have this query, which works well, for an OSM database:

SELECT row_to_json(fc)
FROM (
  SELECT ST_AsGeoJSON(ST_Transform(ST_Intersection(lg.geometry, bbox.geom), 4326))::json As geometry,
  (
    SELECT row_to_json(t)
    FROM (
      SELECT name, type, admin_level
    ) t
  )
As properties
FROM
  osm_admin As lg,
  (SELECT ST_Transform(ST_MakeEnvelope(10.735445022583008, 59.913428634991995,
      10.776472091674805, 59.933989460416115, 4326),3857) as geom) as bbox
    WHERE ST_Intersects(lg.geometry, bbox.geom)
  ) fc;

I'd like to filter this generic query a bit, depending on the zoom level. E.g. I'd like to add filter for zoom level WHERE admin_level < 4. I also need to add other filters for other queries, based on this template (where admin_level is a property of the osm_admin table).

Where can I insert such a WHERE clause?


texts: Enabling layout measuring tool / ruler in ArcGIS Desktop?

texts:
According to ArcMap help, I'm supposed to be able to open the Measure Layout Tool in the Production Cartography toolbar:

Customize > Toolbars > Production Cartography.

described here:

https://desktop.arcgis.com/en/arcmap/10.3/guide-books/extensions/production-mapping/the-measure-layout-tool.htm

However, my Arc version (10.4) does not have this toolbar. I simply can't find a layout measure tool, which is completely ridiculous. The measure tool that I can find only works for map elements, but not for layout elements such as text boxes, data frames, etc.

If you can't measure layout elements how do you create decently laid out maps?

I'm sure I'm missing some simple step or connection. . . though the difficulty I am facing while looking for this information (google searches, tutorial videos, etc.) plus the unavailability of this tool (or any measure / ruler tool in cartographic toolsets) has driven me to ask for advice.


texts: Preventing symbols from different layers from interfering with each other in ArcGIS online

texts:
I am tackling a project that contains several new elements for me since I am a high school social studies teacher and not a full-time GIS professional.

For this particular project I am trying to show the various US interventions during the Cold War around the world.  I have created individual points for each intervention.  These points are in a shapefile that contains a list of 49 different interventions.  The fields include the country of the intervention, the year, a short description and whether or not a certain type of action occurred.

Example:  Guatemala, 1954, (description) and the following values for the six types of action possible: Command=1, Troops=0, Bombing=1, Naval=0, Nuke=1, Jets=0.  What this means is that in 1954 we used Command Operations, Bombing and Nuclear Threats to intervene in Guatemalan affairs.  On this particular point then I would like to see the three distinct symbols for those types of action.  Currently I can achieve this by creating 6 layers but the symbols are all stacked on top of each other.

This has all been done in Desktop 10.1 Advanced.  I plan to create a feature service on our AGOL account once I have solved this problem.

Thanks for any advice.


texts: Error 'DBMS Table not found' when sharing some shapefiles as Runtime Content

texts:
I have a collection of shapefiles that I want to "Share as > ArcGIS Runtime Content". In this collection there are a few shapefiles that cause an error. If I access this function through the menu a "General function failure" is returned. Capturing the error in an ArcPy script is more informative:

... could not be compiled as a Runtime Geodatabase Package. ERROR 999999: Error executing function. Operation Failed
DBMS table not found
Failed to execute (CreateRuntimeContent)

I have looked up other posts involving "DBMS table not found" What does 'DBMS table not found' message mean? and DBMS Table Not Found with SQLite and ArcGIS 10.3) and it seems this issue is associated with vectors shaped in database formats. Since my vectors are all stored in ESRI shapefiles, I am not sure why this error occurs.
The shapefiles seem to be intact when I load them in a map.


texts:
I am currently working on Eastern-European microdistricts, and I want to obtain the addresses that can be found in these microdistricts.

Is it possible to use ArcGIS to solve this issue? I thought about using an OSM Basemap and draw polygons where these microdistricts locate, and after that I can extract the addresses from it, but I am not so familiar in this topic.


texts: Translation Vector effect in the 2D transformation of GPS Network into a local system

texts:
For a project of Bridge monitoring using GPS, (X,Y) GPS measurements have taken and they must be transformed into a local system (x,y) to determine the deformation as shown in the sample picture:

For that, many steps of calculations using some formulas in Matlab are done and summarized as follows:

Determination of the parameters of the transformation.
Transformation calculations of (X,Y) from GPS into a local system (x,y).

After getting the new coordinates (x,y), we performed new calculations to get (x,y) again but with considering the translation vector (provided by Land Authority).

The results we get for both cases are as follows:

Case 1: without translation vector:

Scale factor = (99.8 PPM), Rotation angle = (6.2 ARCSEC)
x(m) = .......,340   , y(m) = .....,992

Case 1: with translation vector:

 Scale factor = (3.2 PPM), Rotation angle = (6.18 ARCSEC)
 x(m) = .......,343   , y(m) = .....,993

From the results, we can see that both coordinates are very similar (differences in mm), but the scale factor is highly affected.

My Question:
Does the translation vector does not has any effect on the coordinates in this type of transformation? Just effect the transformation parameters?

Or it has an effect on the coordinates in other types of transformation?


texts: Building GEE timeseries for multiple locations

texts:
I want to build a timeseries chart for multiple points.
Is it possible to use ui.Chart.image.series and have different lines for different points?
They are not regions, so I don't think .seriesByRegion is appropriate.
When I feed in a feature collection of points, it seems to pull the MODIS data and stretch it along as if it were one location, not parse it out by location like I want.
Code below.
MODIS data for my time of interest
var terraCollection = ee.ImageCollection("MODIS/006/MOD10A1").select(['NDSI_Snow_Cover']).filter(ee.Filter.date('2019-11-01', '2020-05-01'))

Feature Collection of Points that I want to build a timeseries for
var points = ee.Geometry.MultiPoint(
    [[13.427963679755548, 60.28846585303612],
     [14.021225398505548, 60.87127944326994],
     [14.087143367255548, 60.08089287233658],
     [13.120346492255548, 60.02605006413671]]);

Current code for timeseries plot that seems to be merging all the points together:
var testChart = function(dataset) {return ui.Chart.image.series({
  imageCollection: dataset,
  region: points,
  }).setOptions({
  interpolateNulls: true,
  lineWidth: 1,
  pointSize: 3,
  title: 'NDSI over Time at a Single Location',
  vAxis: {title: 'NDSI'},
  hAxis: {title: 'Date', format: 'YYYY-MMM', gridlines: {count: 12}},
  series: {
    0: {color: '#FF0000'}, // Point 1
    1: {color: '#FFFF00'}, // Point 2
    2: {color: '#00FF00'}, // Point 3
    3: {color: '#0000FF'}, // Point 4
  }
})}

print(testChart(terraCollection))


texts: Geometry Generator using polygon boundary expression only returning Exterior Ring QGIS

texts:
I am using a custom QGIS expression function that creates a zigzag line.
See the function here: Wavy Edge on polygon in QGIS using geometry generators
The function needs a line string and not a polygon so I use the boundary expression like this:

This works except if it has interior rings.
If the polygons have interior rings I get this error

Eval Error: MultiLineString geometry cannot be converted to a
polyline. Only single line or curve types are permitted.

And it only shows on polygons without rings.

If use this expression instead:
collect_geometries(make_zigzagline(exterior_ring( $geometry),1.5,3), make_zigzagline(interior_ring_n( $geometry,1),1.5,3))

It only shows those polygons with at least one ring and doesn't show the second ring if it exists.

Is there some way I can have it figure out how many rings it has an use that?
Or is there something I can change in this function I can change to have it accept  MultiLineString geometry?
from qgis.core import qgsfunction,QgsExpressionContextUtils,QgsExpression,QgsProject,QgsPoint,QgsGeometry
@qgsfunction(args='auto', group='Custom', usesGeometry=False, referencedColumns=[])
def make_zigzagline(geom,dist,offset,feature,parent):
    """

converts a linestring to a zig-zag line

    make_zigzagline(geometry,distance(s),offset)

geometry	linestring geometry

distance(s)	linear point distances (single number or a string of comma separated numbers)

offset	perpendicular offset

    Examples:

make_zigzagline($geometry,'15,30',15) -> zig-zag line

make_zigzagline($geometry,15,15) -> zig-zag line

    Use smooth function to create wavelines:

    Example:

smooth(make_zigzagline($geometry,'15,30',15),3)

    """

    if not type(dist) is str:
        dist = str(dist)

    dist = [float(n) for n in dist.split(',')]
    l = geom.length()
    dist_sum = 0
    distances = []
    while dist_sum + round(sum(dist),2) < l:
        for d in dist:
            dist_sum += d
            distances.append(dist_sum)

    # interpolate points on linestring
    points2d = [(lambda g: (g.x(), g.y()))(geom.interpolate(d).asPoint()) for d in distances]
    vertices = geom.asPolyline()
    start = (vertices[0].x(),vertices[0].y())
    end = (vertices[-1].x(),vertices[-1].y())

    points2d.insert(0,start) # prepend start point
    points = [QgsPoint(start[0],start[1])]
    i = 0
    n = 0
    b = -90
    for point in points2d[1:]:
        pt1 = QgsPoint(points2d[i][0],points2d[i][1])
        pt2 = QgsPoint(point[0],point[1])
        a = pt1.azimuth(pt2) + b
        pt = pt2.project(offset, a)
        points.append(pt)
        i += 1
        n += 1
        if n == len(dist):
            n = 0
            b = -b

    points.append(QgsPoint(end[0],end[1])) # append end point
    return QgsGeometry.fromPolyline(points)


texts: How can I cluster a small area (say a circle of 20 meters diameter) in OSM to check the different node Ids in the area

texts: I have a node ID(Say Z) which will be in the road/way. How do I find the node IDs around this point Z in certain radius.


texts: What is the SQL in pgAdmin for PostgGIS Shapefile loader GUI Import and Export

texts: I have been using PostGIS Shapefile loader to import and export shapefile, but now I want to ask is there any SQL alternative to this GUI import and export.


texts: EsriLeaflet - Problem with special german projection

texts:
I got a service with EPSG:3068 (http://epsg.io/3068-15949)

My test code to diplay this service on an Leaflet Map with Esri-Leaflet is as following:

<!DOCTYPE HTML>
<html>
    <head>
        <title>Test</title>
        <meta charset="utf-8" />
        <meta name="viewport" content="width=device-width, initial-scale=1" />
        <!--[if lte IE 8]><script src="assets/js/ie/html5shiv.js"></script><![endif]-->
        <link rel="stylesheet" href="assets/css/main.css" />
        <!--[if lte IE 9]><link rel="stylesheet" href="assets/css/ie9.css" /><![endif]-->

            <!-- Load Leaflet from CDN-->
      <link rel="stylesheet" href="//cdn.jsdelivr.net/leaflet/0.7.3/leaflet.css" />
      <script src="//cdn.jsdelivr.net/leaflet/0.7.3/leaflet.js"></script>

      <!-- Load Esri Leaflet from CDN -->
      <script src="//cdn.jsdelivr.net/leaflet.esri/1.0.0/esri-leaflet.js"></script>

      <style>
        body { margin:0; padding:0; }
        #map { position: absolute; top:0; bottom:0; right:0; left:0; }
      </style>
      <script src="//cdnjs.cloudflare.com/ajax/libs/proj4js/2.0.0/proj4.js"></script>
      <script src="https://cdn.rawgit.com/kartena/Proj4Leaflet/0.7.0/src/proj4leaflet.js"></script>
    </head>
    <body>
        <div id="map" />
        <script>
        var crs = new L.Proj.CRS("EPSG:3068", "+proj=cass +lat_0=52.41864827777778 +lon_0=13.62720366666667 +x_0=40000 +y_0=10000 +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 +units=m +no_defs", {
        origin: [-5469675.610177,5201612.460189],
        resolutions: [
          2.6458386250105836,
          1.3229193125052918,
          0.6614596562526459,
          0.26458386250105836,
          0.13229193125052918,
          0.06614596562526459,
          0.03968757937515875,
          0.013229193125052918
    ]
  });
  var map = L.map('map', {
    crs: crs
  }).setView([52.52261949895406, 13.37957410812372], 5);

  L.esri.tiledMapLayer({
    url: 'http://info.baral.de/arcgis/rest/services/utjsc_de/Grundkarte/MapServer',
    maxZoom: 7,
    minZoom: 0
  }).addTo(map);
</script>

    </body>
</html>

Here is the code on jsfiddle:

Update: New Demo Site: http://s559860669.online.de/EsriLeaflet/index_beta.html

Nothing happens... expect the JS-Console Warning:
L.esri.TiledMapLayer is using a non-mercator spatial reference. Support may be available through Proj4Leaflet


texts:
I'm trying to create my first script with qgis.
I began with graphical modeler but some bugs block me to do calculation on variables.
So I decide to write it myself.

I wrote that for the moment :

##Zonage_alamano=name
##point=vector
##sortie=output vector
##angle=number
##longueur=number
##largeur=number

demilargeur = (largeur)/2

angleinverse = -(angle)

reproj_1 = processing.runalg('qgis:reprojectlayer', point, 'EPSG:2154', None, progress=None)

enveloppe = processing.runalg('qgis:concavehull', reproj_1['OUTPUT'], 0.1, False, False, None, progress=None)

rotate1 = processing.runalg('grass7:v.transform', enveloppe['OUTPUT'],0.0,0.0,0.0,1.0,1.0,1.0, angle, None,-1.0,0.0001,0,None, progress=None)

buffer = processing.runalg('qgis:fixeddistancebuffer', rotate1['OUTPUT'], demilargeur, 12, False, sortie)

But I have this message :

  2017-09-05T16:42:11    1   Erreur : valeur du paramètre faux : None

  2017-09-05T16:42:11    2   Uncaught error while executing algorithm

  Traceback (most recent call last): File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing/core/GeoAlgorithm.py", line 203, in execute
                  self.processAlgorithm(progress)

  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing/script/ScriptAlgorithm.py", line 378, in processAlgorithm
                  exec((script), ns)

  File "", line 21, in
              TypeError: 'NoneType' object has no attribute 'getitem'

I think issue come from v.transform treatment, because I made a test without it and it run well. I searched for information about writting v.transform in qgis script but I didn't find nothing...

Somebody have an idea to help me ?


texts: Removing unwanted Linestrings from Multilinestring in Postgis

texts:
I have a Multilinestring layer and I have noticed that in some of my geometries there are unwanted Linestrings. They are composed of two points (start and end) but both are the exact same points. That is wreaking havoc when i try to read this layer from SQLServer.

Example of one geometry:

   "MULTILINESTRING((-3.00712325415935 43.3437538570123,-3.00711888681593 43.3437506773148),
   (-3.00693096810649 43.34361446115,-3.00692660078307 43.3436112814451),
   (-3.00692660078307 43.3436112814451,-3.00692660078307 43.3436112814451),
   (-3.00692660078307 43.3436112814451,-3.0068846515836 43.3436419720614),
   (-3.00692660078307 43.3436112814451,-3.0069222334601 43.34360810174),
   (-3.00691127288165 43.3436001834041,-3.00690690556026 43.3435970036984))"

Notice third Linestring that has the same start and end points. How could I remove these unwanted Linestrings from my geometries?

One approach could be to remove the linestrings which st_length=0 but I don't know how to build the query. I suppose I need to use st_linemerge first to split the multi geometry into single parts?

I need some expert help on this...


texts: QGIS incoming interaction from web application (web browser)

texts:
Is there a way to instruct QGIS to open and center WFS layer on given feature?

Can we use a link with custom protocol:

Open feature in Qgis

Is there any support in QGIS to handle such or similar links? Will it open new layer every time, or it will be smart enough to center already opened layer  (opening many features one after another) ?


texts:
I am working with QGIS 3.4 and I have written the next code in the Python editor. This aims to create a layer without geometry and save it as CSV file. This works fine but my problem concerns the two last lines: The layer displayed to the layer panel are not permanent despite the CSV file.

from qgis.core import *
from qgis.gui import *
from qgis.PyQt.QtCore import *

root = QgsProject.instance().layerTreeRoot()

# create memory layer
layer = QgsVectorLayer("Point", "tp2", "memory")
pr = layer.dataProvider()

# Enter editing mode
layer.startEditing()

# add fields
pr.addAttributes( [ QgsField("name", QVariant.String),
                QgsField("age",  QVariant.Int),
                QgsField("size", QVariant.Double) ] )

# Commit changes
layer.commitChanges()

# add a feature
feat = QgsFeature(layer.fields())
#feat.setGeometry( QgsGeometry.fromPointXY(QgsPointXY(10,10)))

feat.setAttribute('name', 'john')
feat.setAttribute('age', 10)
feat.setAttribute('size', 180)

pr.addFeatures( [ feat ] )

# Commit changes
layer.commitChanges()

# get absolute file path
projectPath = QDir(QgsProject.instance().absolutePath())
file =  QFileInfo(projectPath,"output.csv")

error = QgsVectorFileWriter.writeAsVectorFormat( layer,file.absoluteFilePath(), "utf-8",driverName="CSV")

if error[0] == QgsVectorFileWriter.NoError:
    print("success! writing new memory layer")
else:
    print("faillure! writing new memory layer,", error)

QgsProject.instance().addMapLayer(layer,False)
root.addChildNode(QgsLayerTreeLayer(layer))

Edit: erratum in the code


texts:
So I've been trying to make a master cost raster (taking into account the slope and landcover data). However I can't make the r.mapcalculator work...
 I tried with r.mapcalc but the interface looks so different from the tutorial I'm using (http://gracilis.carleton.ca/CUOSGwiki/index.php/Evaluating_Landscape_Permeability_in_Quantum) and I ended up not being able to perform the analysis.

I also can't use any of the functions from the "GRASS commands" tab (even though the tab GRASS GIS 7 works fine. I tried installing QGIS over again, taking a look at the directory GRASS is installed and nothing seems to work.


texts: How to determine optimal windmill placement to maximize production?

texts: How can I add automatically the maximum number of points (dots) with 6m diameter and 100m distance each point from other into an area. The dots will present the windmill base (6m) and I want the maximum number of points with 100m distance (Ratio) each in order to reach the highest level of production.


texts: QGIS 3.4.2 - timemanager plugin - NetCDF - attribute

texts:
When updating QGIS from version 2.18 to 3.4.2
my NetCDF files are no longer loading in the timemanager plugin

My time variable in the NetCDF file has the following attribute

double time(time) ;
    time:long_name = "time" ;
    time:standard_name = "time" ;
    time:units = "seconds since 1970-1-1" ;

I have the following message on the log message

INFO    Could not guess timestamp in raster filename. Cause 'NoneType' object has no attribute 'start'
WARNING    Could not import netcdftime. Using fallback computation

I try adding a start attribute without success.
Where could I find information on the timemanager input requirement for NetCDF file?


texts: QGIS get number of unique values within geometry of a different layer

texts:
I have a grid-layer (polygon) and a point-layer. The points that lay within a grid-cell have the same unique id like the grid cell. they also have an attribute "year".
What I want is to create a column in the table of the grid layer that is populated with the number of unique values from the attribute "year" from the point-layer (only of points that lie within the respective cell). So in the screenshot-example, for the upper cell that would be "6" (2014,2015,2016,2017,2018,2019), for the lower it would be just "1" (2 x 2017, so only one unique value), and so on.

In the field calculator I can only find within_geom_count and within_geom_sum, but nothing like within_geom_unique or similar.

I also have found a workaround using the "group stats" add on, where I can count the unique values, and then join the produced csv with the grid layer,  but I'm pretty sure that there is some easy way to do it directly in the field calculator without these extra steps.


texts: How to group by / Order by with logical sequence

texts:
I have a problem in SQL (QGIS). I have a shape with plots. I would like to obtain a certain organisation.
Here a schema of what I have and what I want :


texts: Counting point objects near each feature in polyline layer of ArcMap?

texts:
I am having a little trouble finding a strategy for this. I have a polyline layer with Fiber segments, each have a FiberID. They are not unique, sometimes there will be two segments with the same fiberID.

I also have another point layer that contains the Splice Points. I was wondering if there was a strategy that I could get a count of the Splice points that were within a certain distance, lets say 50 feet of each record in the Fiber Polyline layer. So the result I'd be looking for is a table that shows the number of splices (points) within 50 feet of each record in the polyline table.


texts: Looping through subsets and spatial joining to single output?

texts:
I have two feature class files, point file and polyline file and would like to spatially join pairs by looping/iterating based on a conditional IF statement.   The resulting pairs would then be loaded into the ArcMap spatial join tool (v 10.2.2).  Below is my attempt to block out the logic.

IF Point File Field State Code AND Polyline File Field State Code = X

THEN
Create a temporary subset feature class for each (point and line file) and then load into the ArcMap spatial join tool (Target feature points Join features Lines)

OUTPUT FEATURES CLASS to a new feature class with the spatial join  (that allow all subsets to be inserted).

LOOP (ITERATE)
For the next State Code value until all codes have been matched and the all spatial joins of the subsets are completed.

The idea is to not have to load 50 pairs of point and lines files to complete a spatial join. Spatially joining the whole US set of points and lines at once takes a very large run time. In essence you’ re giving ArcMap a much smaller task at a time only having to spatially join a state of point to a state of lines at a time.


texts: Applying Map Algebra to List of Rasters and Saving New File in New Directory?

texts:
I am trying to multiply each raster by coefficient using Map Algebra in Python. I receive the following error for the code below:

  calcRaster.save(calcRasterName) AttributeError: 'unicode' object has
  no attribute 'save'

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

env.workspace = r"C:\SomeFolder"
outFolder = r"C:\SomeFolder"

arcpy.CheckOutExtension("Spatial")

rasters = arcpy.ListRasters()
for inRaster in rasters:
    calcRaster = inRaster * 720
    calcRasterName = outFolder + "/" + calcRaster
    calcRaster.save(calcRasterName)

Based on some searching,  I redefined the calcRaster variable to:

  calcRaster = (Raster(inRaster) * 720)

However, I receive some comments and a new error:

  File "C:\Users\rrotz\Desktop\Pycharm\RasterCalculator.py", line 13, in

      calcRasterName = outFolder + "/" + calcRaster

  File "C:\Program Files
  (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\sa\Functions.py", line 4329, in
  Plus
      in_raster_or_constant2)

  File "C:\Program Files
  (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\sa\Utils.py", line 53, in swapper
      result = wrapper(*args, **kwargs)

  File "C:\Program Files
  (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\sa\Functions.py", line 4326, in
  Wrapper
      ["Plus", in_raster_or_constant1, in_raster_or_constant2])

  RuntimeError: ERROR 000732: Input Raster: Dataset S:\Group
  Projects\HDF5\2 Tiff Data\8 Monthly Scaled/ does not exist or is not
  supported

The last RunTimeError is regarding the outfolder path that I designated for the final calculated rasters.


texts: Stop click event on input field to be propagated to map layer

texts: I have a simple html text input field within a panel on top of an OpenLayers Map. Everytime I click in the input field, the click event gets propageted to the map. I want the click event not to be propagated to the map when it is within the input field. How can I achieve this?


texts: How to get fields of subLayer of MapImageLayer?

texts:
I have MapImageLayer which has sublayers.
I need to show a popup of sublayer feature and popup should contain all attributes(default popup content) with a custom title. Before 4.12 version we can use below solution:
layer.popupTemplate = {
        title: "SS Gravity Main: {FACILITYID}",
        content: "{*}"
 } as any;

But I am using ArcGIS JS API 4.15 version and now I can't use "{*}" in content as stated in the api-reference doc.
https://developers.arcgis.com/javascript/latest/api-reference/esri-PopupTemplate.html#content
It has 2 solutions:
1 set the Popup's defaultPopupTemplateEnabled property to true - MapImageLayer doesn't support this as stated in api-reference doc.
https://developers.arcgis.com/javascript/latest/api-reference/esri-widgets-Popup.html#defaultPopupTemplateEnabled
2 Set the fieldInfos - SubLayers have a long list of fields. So hard coding these fields doesn't seem a good solution.
The other easy solution is loop through the layer fields to create fieldInfo array:
layer.fields.map(k => {
   return {
     fieldName: k,
     label: k
   };
}

But the issue is I am not getting the sub layer's fields. SubLayer's fields array is null. MapImageLayer has a large no. of sublayers so I can't set outfields as ["*"] for all. PopupTemplate also has outFields property but it is on ArcGIS API for JavaScript 4.9 as stated in api-reference doc:
https://developers.arcgis.com/javascript/latest/api-reference/esri-PopupTemplate.html#outFields
Please suggest any solution to show sublayer's popup with default content or How I can get fields of sublayer.


texts: Multiple styles for the same shapefile in QGIS

texts:
I have a shapefile with the number of inhabitants per area and I should represent both the number of inhabitants and the density. I use to style both with categories.

How do you work in this case? Do you change each times the style of the layer or do you duplicate the layer to have the individual representations for easy check or something else?


texts: I have downloaded Bathymetry data from GEBCO. Its in NetCDF format. I want to create contour from it. I tried so many tools to convert the file into raster but couldn't do so. How can I create raster from netCDF?


texts:
I have several plugins I am using that involve rubberbands and their equivalent points (QgisVertexMarkers).  My problem looks to come from removing them from the display.  This process worked fine in QGIS 3.8, But now QGIS crashes with version 3.10.

The code I'm using is shown here.

    vertex_items = [ i for i in self.iface.mapCanvas().scene().items() if issubclass(type(i), QgsVertexMarker)]
    for ver in vertex_items:
        if ver.data(1)=="COGOStartPoint":
            self.iface.mapCanvas().scene().removeItem(ver)

    m = QgsVertexMarker(self.iface.mapCanvas())
    red = QColor(255, 0, 0,255)
    m.setIconType(4)  #1=cross,  2=x,  3=box,  4=circle
    m.setIconSize(8)
    m.setPenWidth(2)
    m.setCenter(QgsPointXY(pt))
    m.setColor(red)
    m.setData(1,"COGOStartPoint")
    m.show()

It seems to be the first few lines of code that causes the freeze.  Trying to remove the points no longer works, drawing the points is still fine.  I found this technique from PyGIS remove QgisVertexMarker from scene
.

What has changed in version 3.10 from 3.8 that would cause this behavior?

Is this still the best way (the only way) to remove rubberband type graphics?


texts: SLD Polygon Styling: Add Point feature to Northern Edge?

texts:
I am trying to create an SLD to style a polygon; the goal is a simple block fill along with a single point image added to the most northern edge of the polygon.  Something like this:

This example SLD code will add a symbol to the start point of a polygon:

<PolygonSymbolizer>
<Fill>
    <CssParameter name="fill">#000000</CssParameter>
</Fill>
</PolygonSymbolizer>
<PointSymbolizer>
    <Geometry>
        <ogc:Function name="startPoint">
            <ogc:PropertyName>the_geom</ogc:PropertyName>
        </ogc:Function>
     </Geometry>
    <Graphic>
        <Mark>
            <WellKnownName>square</WellKnownName>
        <Fill>
            <CssParameter name="fill">#00FFFF</CssParameter>
        </Fill>
        </Mark>
        <Size>5</Size>
    </Graphic>
</PointSymbolizer>

However I can't work out how to show the point symbol to the most northern edge of the polygon?

One more thing I tried: Using QGIS I can create styles and export to SLD, and using the geometry generator as follows:
     closest_point( $geometry, make_point(  x($geometry), y( $geometry)+20000))
This works.  However when exporting this to SLD it fails due to 'not implemented yet'.


texts:
I usually use ArcMap to display photo attachments in a point layer. If I click the point with the HTML Popup tool, it displays a small photo corresponding to the point and I can click on it to open it in an external view where I can zoom in and out.

I can't find the equivalent in ArcGIS Pro. I can see the attachment using the Explore tool (in the Map tab), but then I can't open it like in ArcMap. I am stuck with the picture the size of a thumbnail.

How do I open the attached photo in an external viewer from the popup windows in ArcGIS Pro?


texts: Grouping and exporting a set number of features in ArcGIS

texts:
I have an areal feature class that has approximately 4,700 features. I need to divide this feature class into multiple feature classes that contain EXACTLY 60 features in each. Location of features have no matter. Yes, that will create nearly 80 feature classes with SplitByAttributes, that is what I need.

There are no attributes to sort with. Grouping Analysis can force creation of 80 groups with attributes, but creates too many features in some files, it is not uniform and not able to force the number of features (that I have knowledge of).

Anyone have any idea how to export a specific number of features to new feature classes?


texts:
I've downloaded OSM file and extracted it to PostgreSQL. Now I have 4 layers:

planet_osm_line
planet_osm_point
planet_osm_polygon
planet_osm_roads

I added these layers to GeoServer and set SLD file for each of them but after rendering road arrow is not shown.
I need something like this:

What should I do to show these arrows?


texts: About the method of resize raster resolution

texts: I opened a raster file as a layer, and click the right mouse button and select "Save as..." to change the resolution. I got the new raster, but I want to know what is the method (average, aggregation...) of this mode?


texts: I simply cannot find out the following: I want to repeat a calculation in GRASS GIS several times. Is there a command for r.mapcalc to repeat a calculation x times (similar to this possibility in r.mfilter, for example)....besides listing that step x times in s shell scripts (I use shell).


texts: QGIS raster area by clipped mask and 2 layers

texts:
I have two layers, one is a tif file with inclination of the soil of a whole county (which I got from the SRTM-relief) and the other file is a shape file which contains all the properties in this same county.
The columns of this .shp file contains area and id.

What I need is to add 3 columns to the attribute table of that shape file with the count of the three raster types by id presented in the tif file. I already managed to count the 3 types of raster for the whole shape.

I uploaded Pictures & Videos here:

https://1drv.ms/u/s!Asdg-BqIM7nshIcXVbPRtHDMyv1mBQ?e=dEhWRu


texts: Ungroup ArcMap Graphic Elements in Layout using ArcPy

texts:
I want to alter text elements in a layout. This text elements were grouped. Grouped elements behave like GraphicElement and it is possible to query the property isGroup. I can not alter these text elements before ungrouping.

Is it possible to ungroup these GraphicElements using arcpy?

EDIT:

My reasoning for this question is false. As Alex Tereshenkov points out: you can alter grouped text elements.


texts: Point plus distance along great circle using RGeo or WKT

texts:
I'm racking my brains (and google) trying to figure out how to get RGeo to do the following geospatial calculation:

Imagine I start off at some arbitrary geographic point (let's say (45, 45) which happens to be in Russia), then I point 30 degrees east of true north, and fly for 5000km, and then want to know my position.

How can I do this with RGeo? Also, is there a way of neatly expressing this in WKT assuming a geographic coordinate system? For instance "POINT(45, 45) + POLAR_VECTOR_THIS_CONSTRUCTOR_IS_MADE_UP(30, 5000000)"

Also, any good references, tutorials, and docs for this sort of thing would be much appreciated.


texts: Connection error with QGIS 3.0 Metasearch plugin

texts:
I've added "Data.gov" as a service in the Metasearch plugin (QGIS 3.0 32-bit). The associated URL is https://catalog.data.gov/csw-all.  When I click "Service info" I receive a response describing the service.  However, when I attempt conduct a search (e.g., "roads") using this service I receive the following error:

             Connection error: mismatched tag: line 4, column 2

Any thoughts on what the problem may be?


texts: Interpolating xyz data in ArcGIS 3D Analyst?

texts:
I'm using ArcScene 10.3 to view temperature data as seen below:

This data come from a CTD, an instrument that measures the conductivity (salinity) temperature and depth of the ocean.  It is lowered to the seafloor making recordings every few seconds.

What I ideally would like to do is interpolate the temperature across the entire area.  Then I could take a vertical slice anywhere and see the expected temperatures.  ArcScene seems great for visualizing the 3D data but not so great for doing any kind of interpolations.

My best solution so far is a total hack.  I choose several of the CDT casts that were somewhat in a line and measured the distance between the first and each subsequent CTD cast.  I then added a new field to the data table called distance with this info.  I used display xy data in ArcMap setting the distance as x and the depth as y and changed the symbol to show the temperature gradient.

From there I interpolated using ordinary kriging and then exported the interpolation as a tif file which I overlayed on the arcScene map.  This solution doesn't seem very rigorous because of the way I calculated distance and the fact that I just placed an image of the interpolation about where it belonged on the ArcScene map. How else can I interpolate this data so I can view in 3D?

*note the first image is from a different area than the last two images


texts: Legend ArcGIS Javascript does not recognize Legend from ArcGIS for server 10.2.2

texts:
I have created maps in ArcGIS for desktop and published them to ArcGIS for server 10.2.2. Because of performance issues I have created a group layer of the same map with different levels of details. They all have the exact same symbology. I add the group layer in my script as one layer. I am attempting to create a legend with it but I am not succeeding. I do get a box with the note "No Legend". This indicates to me it does not recognize my legend from my ArcGIS Server. Could anybody help me?

This question is a follow up of a previous question of me that seems to be fully answered since I am getting the legend box: Add a legend in an ArcGIS Javascript API


texts: Leaflet: get the color of a specified pixel on the map

texts: I need to get the color of a specified pixel on the map. How can I do it using the current API of Leaflet?


texts: How to get features in Openlayers vector layer

texts: What is the exact API to get all features in a vector layer in OpenLayers? Suppose I have the "map" object, within which is a "Vector Layer": how can I obtain features in the layer?


texts: Can't view RGB image exported to Google Drive from Google Earth Engine

texts:
I am trying to export a clipped satellite image from Google Earth Engine's Code Editor. I am able to see the image correctly within code editor, but am unable to open the exported image in Google Drive or on my computer if I download it.

My code is below, does anyone know what the problem might be?

/**
 * Function to mask clouds using the Sentinel-2 QA band
 * @param {ee.Image} image Sentinel-2 image
 * @return {ee.Image} cloud masked Sentinel-2 image
 */
function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

  return image.updateMask(mask).divide(10000);
}

var geom = ee.Geometry.Polygon([[33.8777, -13.4055],
                                [33.8777, -13.3157],
                                [33.9701, -13.3157],
                                [33.9701, -13.4055]])

var collection = (ee.ImageCollection("COPERNICUS/S2")
                  .select(['B4', 'B3', 'B2', 'QA60'])
                  .filterBounds(geom)
                  .filterDate('2016-01-01', '2016-06-30')
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                  .map(maskS2clouds)
                 )

var image = collection.sort('system:index', false).mosaic()
var image = image.clip(geom)

// Export the image, specifying scale and region.
Export.image.toDrive({
  scale: 10,
  image: image,
  description: 'example_image',
  fileFormat: 'GeoTIFF',
  region: image.geometry().bounds().getInfo()['coordinates']
});


texts: Converting raster from LCCS to NLCD formats

texts:
We would like to convert LCCS format raster files into the US NLCD format. Is there some way for that, especially with QGIS program ?

We need NLCD TIFF files for EPA AERSURFACE program. The https://lcviewer.vito.be/ world cover has only LCCS land cover data.


texts: How to check, if ArcMap custom tool is active? I have MyCustomTool and docWindow. docWindow opens when MyCustomTool is created. When I close docWindow and select another tool (hand - for example), then MyCustomTool again - docWindow doesn't open.


texts:
I am trying to project APJ specific meridian on -180, under WGS84/PDC Mercator always produces gap as attached.

Can you help explain how to manage an overlapped projection?


texts: Labelling only the selected features with custom text using pyqgis

texts:
I want some alternate to this code as this code just labels only the current feature in the loop .

FEATURES = point_layer3.selectedFeatures()
for  f in FEATURES:
    total = f['parkeer_approx']
    alloted = parking_capacity[f.id()-int(1)]
    print(str(f['identificatie']), alloted['parking_id'])
    layer = point_layer3
    layer.setCustomProperty("labeling", "pal")
    layer.setCustomProperty("labeling/enabled", "true")
    layer.setCustomProperty("labeling/fontFamily", "Arial")
    layer.setCustomProperty("labeling/fontSize", "10")
    layer.setCustomProperty("labeling/isExpression", "true")
    layer.setCustomProperty("labeling/fieldName", 'CASE WHEN "identificatie" = \'{}\' THEN "parkeer_approx" ELSE '' END'.format(str(f['identificatie'])))

    layer.setCustomProperty("labeling/textColorG",255)
    layer.setCustomProperty("labeling/placement", "2")
    iface.mapCanvas().refresh()

What I wanted to do is that I want to label all the features if they satisfies the CASE. Suppose there are 8 features selected in the layer so i wanted to label only those 8 features one by one if they satisfies the CASE.


texts: Snapping dangling nodes with prescribed angle using ArcGIS Desktop?

texts:
I have got a streets shape file, to snap the dangling nodes in ARCGIS.

But, I have to snap only some nodes(of feature) with a prescribed angle (green node in the image, and the angle marked in grey). And the other features (marked in red) has to be left out.

Any preset tool or python script available?


texts:
I can't remove a Leaflet layer from my map. I have referenced a lot of other SO threads about this very subject. Even so, removeLayer() refuses to work. addLayer() works just fine.
Threads I've visited (among others):
leaflet-layers-on-different-zoom-levels-how
leaflet-how-do-you-use-removelayer
function-with-map-addlayer-and-map-removelayer-not-working
map.on('zoomend', function() {

    const npsBoundaries = L.esri.featureLayer({
        url: 'https://services1.arcgis.com/fBc8EJBxQRMcHlei/ArcGIS/rest/services/National_Park_Service_Boundaries/FeatureServer/0',
        style: function(feature) {
            return { color: 'gray' };
        }
    });

    if (map.getZoom() > 6.5) {
        map.addLayer(npsBoundaries);      // <-- this works
    } else if (map.getZoom() <= 6.5) {
        map.removeLayer(npsBoundaries).   // <-- why doesn't this work?
    }

If it helps, I am using Leaflet v1.7


texts: Migrating from CARTO Editor TO CARTO Builder and problem with odyssey.js (viz.json)

texts:
In odyssey.js you have:

-vizjson: URL like this "//username.cartodb.com/u/usernanme/api/v2/viz/key/viz.json"

The first thing we want to do is bring the CartoDB Visualization we created into the Odyssey.js application. We're going to use this vizjson option in the config block of your Odyssey.js markdown.

But now we cannot get this vizjson URL from cartDB!

The CARTO Builder does not have the option to share the viz.json URL due to the current version of CARTO.js (v3.15) which is not compatible with the Builder.

see:  https://carto.com/learn/guides/intro/migrating-from-carto-editor-to-carto-builder

MIGRATING FROM CARTO EDITOR TO CARTO BUILDER:

Can I use CARTO Editor and CARTO Builder at the same time?

No, once CARTO Builder is enabled for your account, it is the only available map application.

Can I choose to use CARTO Editor instead?

No, CARTO Editor is currently in the process of being deprecated. All of the Editor features, and additional advanced functionality, is available within CARTO Builder. Once CARTO Builder is enabled for your account, the CARTO Editor is no longer available.

All new accounts have CARTO BUILDER enabled...
If somebody could help me with an old account (for example without datasets/maps), it will be nice.


texts: Specifying a file for VirtualShape in SpatiaLite

texts:
I'm fairly very new to spatialite, and am stuck on something very basic.

I'm creating a virtual table from a shapefile using virtualshape. This works fine when the shapefiles are in the same folder as spatialite.exe.

CREATE VIRTUAL TABLE test USING VirtualShape(test, CP1252, 28356);

I can't figure out how to specify an absolute path ie C:\1. Project Data\1814\test

CREATE VIRTUAL TABLE test USING VirtualShape("C:\1. Project Data\1814\test", CP1252, 28356);

Fails.

What's the syntax that's required here?


texts: I know in QGIS you can subdivide a line, but is there a way you can type in a distance and the line will split according to this, with the last segment being smaller in length? If not I will try to script this in Python.


texts: OpenLayers programmatically select base layer

texts:
I have several base layers added to a map in an OpenLayers based application. There is only one that can be displayed at a time, as they are all base layers.

When the user changes the selected base layer, I store the name of the selected layer as a HTML5 local storage item (an improved cookie). When re-opening the application, I want to programmatically set the same base layer be displayed.

It works fine with layer.visibility = true if isBaseLayer is false, but not in the case when isBaseLayer is set to true.

Is there anything else I could try for this?


texts: Splitting Trimble SSF file for post processing 2 days separately?

texts:
One of my field guys collected positions (updated features in an .imp file that I generated from our billing system) in the same file over multiple days. Usually I "split" the file at the end of each day for him (extract features from the file including updated positions), but I got really busy on another project and forgot this past time. Of course, this would happen over the weekend of the time change.

Now I have an SSF I cannot differential correct because it spans 2 time zones.
I need to split this file into the 2 days it was collected, AND I assume change the date/timezone in the header in order to correct each day separately. How can I achieve this?

If I use SSF Record editor, delete (strike through) the positions from one day and save-as, I assume the post processing will still fail since the collection range in the header is across multiple days including a time change.


texts: Sum multiple rasters using raster calculator, the result raster had a larger cell size?

texts: I have 365 rasters (daily) and therefore I summed them to form an annual raster using raster calculator in Arcgis. The original cell size of daily raster is 5*3.75 (col 72, row 48), while the cell size of the annual raster is 5*5 (col 72, row 36). The cell size became larger automatically. Does anybody know why?


texts: Creating heat map based on attribute from another layer in QGIS

texts:
I have a layer showing different trails in a nature areal. For each trail, there is an attribute with the usage of the trail (no. of people passing a trail counter). I now want to create a new layer as a heat map, based on the attribute. How do I do that?
I only know how to do in the original trail-layer, working with the symbology. But I want the heat map to be showed as a separate layer, since it's easier to work with the original trail-layer. So, the heat map should automatically load the geometry and attributes from the original trail-layer.


texts: Can Geoprocessing Services be consumed within ArcGIS Online?

texts:
I have a geoprocessing package in my ArcGIS Online account for my organization which can be downloaded and used within ArcGIS for Desktop.

Are similar functions available to enable Geoprocessing Services be consumed within the ArcGIS Online environment?


texts: Difference between a spectroradiometer and multispectral sensor

texts:
Let's take as an example of a spectroradiometer MODIS, and as an example of a multispectral sensor the Landsat8 (OLI).
What is the difference between these two sensors?

Should we compare a spectroradiometer instead with a hyperspectral sensor since in the case of MODIS the number of bands is more important (36 bands).


texts: Creating QGIS action which copies features between layers?

texts:
i tried (i never tried before...) to write a simple python code like action in QGIS on my "XXX" vector layer that:

open a layer "YYY"
copy the selected feature on my "XXX"layer
paste the selected feature in the "YYY" layer

Point 1) work ok, but point 2) don't work

Here is my code:

qgis.utils.iface.addVectorLayer("/home/stefano/Documenti/YYY.shp", "YYY", "ogr");
qgis.utils.iface.actionCopyFeatures()


texts: Why Python says shapefile does not exist?

texts:
I'm working on a script that prints the names of hospitals that are within a certain zip-code. When I run my code, though, python keeps giving me an Error message that says that my Hospitals shapefile does not exist. I've tested to see whether or not the shapefile exists using arcpy.Exists and I get the answer "True". I was also wondering if it might be because the file is locked in some way but I've tried deleting my objects and cursors as well as closing arcgis and restarting python and nothing has changed.

Can anyone tell me why Python is not recognizing my shapefile?

Code for reference:

import arcpy
arcpy.env.workspace=r"C:\Users\Cara\Documents\MPSGIS\Programming and Scripting\Lab6\Lab 6 Data"
work=arcpy.env.workspace
zip5=78705
strZip=str(zip5)
arcpy.env.overwriteOutput = True

def searchHospitals(work,zip5):

    fc= "Hospitals.shp"
    #create a list to store values
    hospitals=[]

    whereClause= '"ZIPCODE" LIKE'+"'%s%%'"% strZip
    sCursor= arcpy.SearchCursor(fc, whereClause, "", "", "")

    for row in sCursor:
        rowHospital = row.getValue("NAME")

        #if rowHospital value is not already in hospitals, append.
        if (not(rowHospital in hospitals)):
            hospitals.append(rowHospital)

        #number of hospitals in list
        count= len(hospitals)

        #print results
        print hospitals
        print "Total %s hospitals were found in zip-code %s"%(count, zip5)

    del row, sCursor

Here is the error message:

searchHospitals(work,zip5)

   Traceback (most recent call last):
  File "<pyshell#4>", line 1, in <module>
    searchHospitals(work,zip5)
  File "C:\Users\Cara\Documents\MPSGIS\Programming and Scripting\Python Scripts\cursor_objects.py.py", line 98, in searchHospitals
    sCursor= arcpy.SearchCursor(fc, whereClause, "", "", "")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\__init__.py", line 1167, in SearchCursor
    return gp.searchCursor(dataset, where_clause, spatial_reference, fields, sort_fields)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\geoprocessing\_base.py", line 359, in searchCursor
    self._gp.SearchCursor(*gp_fixargs(args, True)))
IOError: "Hospitals.shp" does not exist


texts: Importing 16-bit data into Google Earth Engine

texts:
I am trying to import 16-bit data into Google Earth Engine.

I get this error.

  No data value 65536.0 cannot be applied to band #0 of type Integer<0,
  65535>.

No Data is 65536

I tell it that 65536 is no data here


texts: Extract data from NetCDF using a shapefile polygon in Python

texts:
I need to subset data from a NetCDF using a specific shapefile.  The data are sea surface temperature and ocean color at 1/4 degree resolution.  I have 4 polygons describing the US. Northeast continental shelf large marine ecosystem and it's sub-components that I need to use to extract the data.  I am working with monthly composite files from 1982-2014, so this data extraction routine needs to be automated.  The files are already subsetted to the approximate working area grid of [35, 45, -80, -60].

Previously, we were converting HDF5 data files to rasters in R and processing them this way, but this method is really inefficient and I am sure there is a better solution in Python using the current NetCDF files.

Thus far I have been using GDAL and Fiona to read in the shapefiles and NetCDF4 to load the data files. I am not sure how to go about subsetting the data.  I found this:

GDAL for Python: extracting subdomains from NetCDF file?

But I don't have the foggiest idea about how to subset a NetCDF file using anything other than a simple bounding box, which these polygons most certainly are not.

Point in polygon routines would probably take an eternity to work, but maybe I could subset the data using a smaller bounding box that is rotated to fit these shapes like this as an initial starting point and then do a point-in-poly search:

Subsetting a curvilinear netCDF file (ROMS model output) using a lon/lat bounding box.

Any ideas?

EDIT 1:

I just came across the OpenClimateGIS package which seems that it may fit the bill perfectly...  I will have a go with this to see if I can get it to work:
http://ncpp.github.io/ocgis/examples.html#advanced-subsetting


texts: Error "Extension is invalid for the output raster format" while batch projecting .sid files

texts:
I'm attempting to batch project .sid rasters using python. I'm able to successfully use the tool in arcmap so I know the software is capable of doing it. Here is my code:

try:
    import arcpy, sys

    InFolder = r"C:\Original"
    OutFolder = r"C:\Reproject"
    OutSR = arcpy.SpatialReference(26918) # NAD83 / UTM zone 18N

    arcpy.env.workspace = InFolder

    for Ras in arcpy.ListRasters():
        arcpy.AddMessage("Projecting " + Ras)
        arcpy.ProjectRaster_management (InFolder + "\\" + Ras, OutFolder + "\\" + Ras,OutSR)
    arcpy.AddMessage("Projecting complete")

except:
    print "Project Raster failed"
    print arcpy.GetMessages()

I'm getting the following: "ERROR 000445: Extension is invalid for the output raster format."

I'm thinking it has to do with the .sid format being turned into a .tiff but I'm unsure about how to get the program to output .tiff format.


texts: Google Earth Engine export multiples images in one task

texts:
I am facing to a problem caused by the number of requests that can be done with the Google Earth Engine Python API is limited to 3000.

Since I have a lot of small disparate images with different regions to export I reach this limit every time. Then to reduce the number of requests done to the Python API I would like to clip these images in one task. The problem is each image corresponds to a different region and following the documentation it seems that it is not possible to get a list of region linked to images clipped.

A pseudo-code based on the existing code would be:

Export.image.toDrive({
  image: [img1, img2, ..., imgn],
  folder: "name",
  description: 'descr',
  scale: 30,
  region: [geometry1, geometry2, ..., geometryn]
});

Is it a way to achieve that?


texts: Handling multiple outputs in QGIS 3.0 Processing plugin?

texts:
I'm porting my QGIS processing plugin algorithm from 2x - 3x.

The algorithm generates and outputs one to many rasters and adds them to the QGIS map canvas. The number of output rasters is unknown until the algorithm is executed.

In my QGIS 2x version, I set the output of the algorithm to be a processing.core.outputs.OutputDirectory,  save all the output rasters in that directory, then add them to the map using the QgsMapLayerRegistry.instance().addMapLayer method.  This is a bit of a kludge as it means the outputs can't automatically be used as inputs to other processing algorithms, but the 2x version was just a quick proof of concept demo.

Other than doing something similar to my 2x version (i.e. set a QgsProcessingParameterFolderDestination/QgsProcessingOutputFolder and add them to the map via QgsProject.instance().addMapLayer), is there a recommended approach for QGIS 3x processing algorithm plugins to handle an unknown number of outputs? Perhaps something analogous to an arcpy.Parameter with direction='output' and multiValue=True properties.

I would like to add the output rasters to the map when the alg is run manually from the toolbox (which I can), but this is not the primary aim of this question. It's more about having the multiple output rasters usable as inputs to other algorithms in a workflow.


texts: Where do I find the legend for netCDF files in QGIS?

texts:
I'm using QGIS to make maps of the maize, rice and wheat cultivated area at global level. I found these great netCDF files: http://www.geog.mcgill.ca/landuse/pub/Data/175crops2000/NetCDF/

Here's an example of the maps that they can generate: Maps

I have loaded the maize, rice and wheat files as layers in QGIS, and get a rainbow of colors. Problem is, I have no idea what they mean. According to the legend on the map above, the color scale indicates "% of total area"; the same datasets were used - I think - for this map, where the legend indicates "Cassava area (Percent cropland)".

Both those legends are bit cryptic (and the colors I get in QGIS are around 40). So I'm trying to locate where in QGIS I can find the "official" legend and descriptions of what the colors indicate. I followed the steps here and the legend I got was just the names of the three layers, with no colors. If I try just one layer, I get the name of the layer and nothing else.


texts: Exporting polygon shapefile with lat/lon points

texts:
By copying rows of the attribute table of a polygon shapefile, I figured out how to extract the points of the given polygon but to my surprise, they are formatted backwards:

POLYGON((
-87.529571533203125 34.55934906005859375,
-87.52953338623046875 34.3049163818359375,
-87.93964385986328125 34.31328964233398438,
-88.17311859130859375 34.32147979736328125,
-88.1404571533203125 34.58187484741210938,
-87.95177459716796875 34.57530975341796875,
-87.52912139892578125 34.56753158569335938,
-87.529571533203125 34.55934906005859375))

Is there a way to export these points in the same way but formatted as lat/lon xx.xx, -xx.xx points?(only US counties are being processed).
If not and post-formatting is the only option, if there is anyone out there who has already done this...how?
If not I guess I will just have to write a formatting script but it will be a pain to execute on all 4703 counties.


texts: Set Browser location dynamically to random location

texts: I am working on a program that is supposed to run on a mobile client and track the location. However, as I am implementing on a non-mobile desktop client, I would like to have something like a chrome plugin that changes my location dynamically within specified time intervals (or shuffles my location somehow else). Is it know if something like this exists?


texts:
I'm working with the javascript library OpenLayers 3.0. I try to have the listing of the layers of my map.
If I do
console.log(map.getLayers());
I've got the array [a] my four layers :

But I don't find in that, the name of my layers.
How can I take the name and print them in an alert ?


texts: Finding good Mbtiles to png tiles converter?

texts:
I'm trying to create a slippy map and I need to obtain png tiles. I have exported a .mbtiles file from TileMill and I tried several python scripts to convert this file into png tiles but I always end up with a messy folder where the tiles are misnamed, and I can't use them into my app.

Does anyone know a good converter which organize and name your tiles properly?


texts:
Here is a pretty pic of the raster I am using - with no effects applied - just the color ramp and histogram adjusted to show out all the relief - albeit in a very 2D way!

I have a raster layer that I want to display with hillshade as a default.  I can do this in code and in effect it is the same as selecting the Hillshade option in ArcMap (see below):

When you do this however - it looks disappointing:

Especially if you zoom into the raster you see all kinds of gridded articfacts - which does not look nice at all!

However if I create a proper Hillshade layer using the raster with Spatial Analyst and overlay it with transparancy of 85% it looks pretty good and very different!:

How can I acheive this effect in ArcMap and override the default functionality of the Hillshade option (I supsect this is impossible, or is it some setting I have missed out that is caused the raster to be displayed poorly).  If it is possible to override or select the correct settings - is it possible to do in ArcObjects.
If option 1 is not possible, can I create my effect with creating a Hillshade layer (without using Spatial Analyst) and effectively fuse it to the orginal raster - i.e. I don't want people to see two layers in the table of contents in ArcMap. I just want one layer but with a default hillshade rendering of my choosing/configuration.  I need to do this programmatically too i.e ArcObjects.

One thing I have learned is that you must have your DEM in the correct Spatial Reference and measurement units or the resulting hillshade looks too black and over exaggerated - check out this video to find out more - its a very good explanation.

So I think I now know what I need to do in regards to producing some high quality hillshading, but I think it's going to take some low level tweaking ArcObjects code.  I am thinking I need to come up with my own graphics rendering and plug this into ArcMap (I have a dll that renders bitmaps of grids in a cool way).  I think a good place to show this would be to override the Draw() event on something like the IRasterRenderer, and insert my own image rendering.  Of course it might get a bit complicated with reprojections etc.

Has anyone done any kind of this low level graphics rendering in ArcMap - if so can you offer any guidance or things I should be aware of?


texts: How to select a shape and round the attribute field?

texts:
# Import arcpy module
import arcpy

if arcpy.Exists("F:\\Atco_gis\\temp_shapefiles\\Temp_LP_Layers\\adverse_skid_area.shp"): arcpy.Delete_management("F:\\Atco_gis\\temp_shapefiles\\Temp_LP_Layers\\adverse_skid_area.shp")
if arcpy.Exists("F:\\Atco_gis\\temp_shapefiles\\Temp_LP_Layers\\comments.shp"): arcpy.Delete_management("F:\\Atco_gis\\temp_shapefiles\\Temp_LP_Layers\\comments.shp")
if arcpy.Exists("F:\\Atco_gis\\temp_shapefiles\\Temp_LP_Layers\\desig_crossing.shp"): arcpy.Delete_management("F:\\Atco_gis\\temp_shapefiles\\Temp_LP_Layers\\desig_crossing.shp")

# Local variables:
adverse_skid_area = "F:\\Atco_gis\\temp_shapefiles\\Temp_LP_Layers_GDB\\Temp_LP_Layers.gdb\\adverse_skid_area"
comments = "F:\\Atco_gis\\temp_shapefiles\\Temp_LP_Layers_GDB\\Temp_LP_Layers.gdb\\comments"
desig_crossing = "F:\\Atco_gis\\temp_shapefiles\\Temp_LP_Layers_GDB\\Temp_LP_Layers.gdb\\desig_crossing"
Temp_LP_Layers = "F:\\Atco_gis\\temp_shapefiles\\Temp_LP_Layers"
Temp_LP_Layers__2_ = Temp_LP_Layers

# Process: Feature Class To Shapefile (multiple)
arcpy.FeatureClassToShapefile_conversion("F:\\Atco_gis\\temp_shapefiles\\Temp_LP_Layers_GDB\\Temp_LP_Layers.gdb\\adverse_skid_area;F:\\Atco_gis\\temp_shapefiles\\Temp_LP_Layers_GDB\\Temp_LP_Layers.gdb\\comments;F:\\Atco_gis\\temp_shapefiles\\Temp_LP_Layers_GDB\\Temp_LP_Layers.gdb\\desig_crossing", Temp_LP_Layers)

arcpy.CalculateField_management("F:\\Atco_gis\\temp_shapefiles\\Temp_LP_Layers\\adverse_skid_area.shp","PERIMETER","!Shape_Area! / 10000","PYTHON")

round(!PERIMETER!, 2)

How do I select the adverse_skid_area.shp perimeter attribute and round the field to 0.00 (two decimal places?)


texts: How to identify train stations from (vague) codes?

texts:
I'm trying to identify and geocode several train stations based on location codes that were supplied to me. Most of them take the form "CCSSS" where "CC" appears to be an abbreviation for a country (such as "GB" or "BE" or "US") and "SSS" appears to be a station code.

However, some of the station codes don't match either an IATA code, a British station code, or an Amtrak station code. Some of them are only three letters long. Some of them contain an "@" symbol. The full list is uploaded here. Interestingly, Google was able to identify "GBQQS" as St. Pancras when I searched for "QQS train station great britain".

How can I reliably identify these stations?

If it's relevant, I'm planning to geocode the stations using the Google Maps API.


texts:
I am trying to rasterize a tif file by doing this: newData <- raster("file.tif"). However, when I run it, it says that it "Cannot open TIFF file due to missing codec". Upon looking online, it seems like the solution to this problem is to either update or install libtiff. I'm on Windows. How do I install libtiff into RStudio?
Things I tried:
I downloaded libtiff-0.4.2.tar, and am trying to install into my RStudio via install.packages("libtiff-0.4.2.tar", repos = NULL, type = "binary"). However it gives me this error
Warning in install.packages :
  error 1 in extracting from zip file
Warning in install.packages :
  cannot open compressed file 'libtiff-0.4.2.tar/DESCRIPTION', probable reason 'No such file or directory'
Error in install.packages : cannot open the connection

I also did install.packages("tiff"), library ("tiff") in RStudio.
I also installed it in Anaconda via "conda install -c anaconda libtiff" but as my Anaconda is synced up with Python, it didn't do anything for R. Is there a way to do the same thing but in Command Prompt/Powershell for R? Sorry, I'm all very new to this and unfortunately have a deadline coming up soon


texts: Problem exporting shp to CAD with ArcMap 10.4 - all features export to layer 0

texts:
I am finding trouble when using the tool to export shapefiles to dwg when using arcmap 10.4.

After using the tool to create CAD fields, I name the fields "Layer" and "RefName" the same way, and set the value of the field LyrOn to 1. As far as I know this is the requirement to create a CAD layer when exporting from arc to dwg.

When running the conversion tool no error messages appear and all seems to go nicely, but when I open the dwg, all the features have exported onto Layer0, and none of the layers specified in the "Layer" and "RefName" fields of the shapefiles appear.

If I use a seedfile where the layers are specified beforehand, the conversion will go alright, but whatever is not specificed in the seedfile, will go automatically to Layer 0.

Before upgrading to version 10.4, whatever wasn't in the seedfile would not export to CAD at all - now it does onto Layer0, so I guess it's progress.
However, clearly something is wrong and I don't quite get what it is.


texts:
I want to remove the direction for latitude (N) and Longitude (W) and keep (°) like 52° and -130° in the map. I tried to find in the draw coordinates and format but I haven't found a way yet. I think I have to write a program in custom but I don't know it.

Edit:
Here the example of the map containing latitude and longitude values with respectively N and W. I just want to remove these two letters.


texts: Run arbitrary Python scripts on GeoServer

texts:
I need to create a Web Processing service on GeoServer (2.11.2) that triggers a back-end Python script. This script should be able to utilize requests and other packages necessary to do the job.

I have got the wps extension and all necessary GeoScript plugins installed with GeoServer.

I followed an example as below to create a custom script that can be called from the Admin Web Console WPS Request Builder.

http://suite.opengeo.org/opengeo-docs/processing/scripting/processcreate.html

import math, sys
from geoserver.wps import process
from geoscript.geom import Point
from geoscript.feature import Feature
from geoscript.layer import Layer

@process(
  title = 'Distance and Bearing',
  description = 'Computes Cartesian distance and bearing from features to an origin.' + str(sys.executable),
  inputs = {
    'origin': (Point, 'Origin from which to calculate distance and bearing.'),
    'features': (Layer, 'Features to which distance and bearing should be calculated.')
  },
  outputs = {
    'result': (Layer, 'Features with calculated distance and bearing attributes.')
  }
)
def run(origin, features):

  for f in features.features():
    p = f.geom.centroid
    d = p.distance(origin)
    b = 90 - math.degrees(math.atan2(p.y - origin.y, p.x - origin.x))

    yield Feature({'point': p, 'distance': d, 'bearing': b})

I tried to use Python's sys library and print out sys.executable as a testing process, but the output was NoneType.

I also tried to import requests but the GeoServer logs showed requests cannot be found.

I am wondering whether arbitrary Python packages e.g. Requests can be used in a script as such supported by GeoServer WPS. Are there any limitations for Python scripts that can be triggered as Web Processing Service by GeoServer?


texts: Get BBOX values from a WMS GetCapabilities response with OpenLayers 2

texts:
I am using the below code to get the BBOX values:

wms = new OpenLayers.Format.WMSCapabilities();
OpenLayers.Request.GET({
url:"http://localhost:8080/geoserver/wms?service=WMS&version=1.1.1&request=GetCapabilities&",
 success: function(e){
  var response = wms.read(e.responseText);
  var capability = response.capability;
  for (var i=0, len=capability.layers.length; i<len; i+=1) {
  var layerObj = capability.layers[i];
  if (layerObj.name === 'cite:'+layers) {
  var bounds=OpenLayers.Bounds.fromArray(layerObj.bbox).transform(new OpenLayers.Projection("EPSG:4326"), new OpenLayers.Projection("EPSG:3857") );
    mapPanel.map.zoomToExtent(bounds);
    break;
  }
  }
  }
});

But I am not getting the exact values, when I debug the code in Firefox it look like this:

I want the BBOX values shown in the above image. Please tell me the code to get that BBOX values.


texts: Catching a system exit caused by ArcPy in PyCharm?

texts:
Using ArcPy in PyCharm, I'm constantly getting Process finished with exit code -1073741819 (0xC0000005) in my RegionGroup loop. I have posted my code in a previous question (Process finished with exit code -1073741819 (0xC0000005) from RegionGroup (ArcGIS 10.7)).
Here I want to know if there is a general possibility to catch a system exit like this (by transforming it into a "real" python error maybe), skip the current iteration and start the next.

I've looked at Catching a system exit with GDAL Python but this seemed to have been due to a GDAL-related bug.

EDIT:
I've tried catching the error using try and except, but so far that has not worked.

try:
    regionGrp = RegionGroup(sinks)
    # or to avoid saving in a tmp file with potential access violations:
    region_save = "C:/.../regionGrp{}.tif".format(glacierID)
    RegionGroup(sinks).save(region_save)
except:
    print "== strange system exit caught =="
    break

The result in the terminal (if the loop breaks) is:

Connecting regions
Compacting labels

Process finished with exit code -1073741819 (0xC0000005)

Normally it should go on like this:

Connecting regions
Compacting labels
Building Attributes
Updating link item

Maybe the error cannot be caught because it happens in the middle of this multi-step RegionGroup algorithm? I'm fairly new to python so I'm not sure if this could be a possibility.


texts: How to calculate unique IDs starting from highest current ID

texts:
I have been using a python script in field calculator in ArcGIS 10.2 to add unique IDs automatically, however there are two things I would like to improve this process:

Each ID also has letters at the beginning to identify what it is (ex. SL is streetlight) Currently I calculate the ID and then go back into field calculator and add 'SL'+(field) to append the text, is there a way I can put this right in my code?
When I add a new feature I need to check what the current highest ID is so it can go from there. This is made more difficult because when I sort by ascending/descending, the field is text and therefore 12 would come before 2. Is there a way to automatically calculate the highest existing ID?

Code in field calculator-

Codeblock:

rec=0
def autoIncrement():
    global rec
    pStart = #last unique ID
    pInterval = 1
    if (rec == 0):
        rec = pStart
    else:
        rec = rec + pInterval
    return rec

Basin:

autoincrement()


texts: importing data (weighted average) from a shapefile from another shapefile

texts:
I am new to GIS and hope this is not too basic. I am using QGIS.

Say I have two shapefiles: SHP1, with attributes allocated to municipalities, and another shapefile I built from a wkt file, SHP2, which is an orthogonal 60 by 60 grid.

I wish to do the following: add an attribute column to SHP which reflects the weighted average of an attribute from the underlying SHP1.


texts: Accuracy assessment in QGIS (Semi-Automatic Classification Plugin) without using random points?

texts:
i am using QGIS with the Semi-Automatic Classification Plugin. I already classified a Sentinel 2a scene by using an supervised classification approach. So I have trained areas all over the scene but I do not want to have an accuracy assessment for the whole scene, but for defined areas in this scene. I want to check every single Pixel in those defined areas (150 - 250 pixel). Does anybody know wether this is possible and if so:

It is possible within the Semi-Automatic Classification Plugin and how?
or
Do I have to use another tool in QGIS (like Grass, Saga etc.) and which one?


texts: Reduce feature count when polygonize in QGIS

texts: I'm performing a polygonize to convert a raster to a vector. Is there an option to reduce the quantity of features that it picks up from the raster file, at the moment there is 12 million features. I have 12 of these to do then merge them which is too large and slow. I'm using QGIS with GRASS.


texts:
I have a point feature class in a geodatabase that uses WGS84 as its SRS. A colleage of mine created points in this feature class using an external database containing lat and long values. However, the lat and long values in the external database were using NAD83 (my fault -- I knew they were, I just forgot). I thought it wouldn't be a problem -- I just created a copy of the geodatabase and changed the coordinate system used by the feature class to NAD83. Should be good, right?

I added this feature class into ArcMap. I saw that the points in the copied feature class (now using NAD83) exactly line up with the original feature class (using WGS84). That doesn't seem right. If two points, having the same lat / long values, but referenced in different coordinate systems are overlaid on the same map, one of them should be projected on the fly to the coordinate system of the data frame, no?

Now, here's the weird part -- if I project the original point feature class (using WGS84) to NAD83, and add it to the map, then the points don't line up. See, this isn't what i expected, either. I expected that they would be projected to the other coordinate system (i.e., their lat / long values would be different by virtue of being projected, but they would be aligned with the original point because they were simply projected to the new coordinate system (then projected on the fly back to the original coordinate system).

Can anyone explain the behaviour that I'm seeing? And, how do I properly ensure that the NAD83 lat and long values collected by my colleague are projected to WGS84 in my target geodatabase properly?


texts: Programmatically converting arbitrary XML data to shapefile?

texts:
Based on the answers below I decided to do some programming. I'll be using GEOtools for this, a java lib: http://www.geotools.org/

I have the following XML:

<shapes xmlns="http://www.meteogroup-maritime.com/spos/GISLayer"
    name="Load line zones" transparency="50" onland="false">

 <shape id="0" description="SUMMER ZONE" color="#FFFFC90E">
   <polygon>
    <location lat="35" lon="-180"/>
       ....
    <location lat="-33" lon="-170"/>
    <location lat="-47" lon="-180"/>
   </polygon>

   <label lat="29" lon="-45" text="SUMMER ZONE" />
   <label lat="-30" lon="-153" text="SUMMER ZONE" />
 </shape>

  ...

</shapes>

I want to turn this data into a shapefile, preferably by using existing software, although I don't mind scripting.

How can I do this?


texts: Color polygons in different colors using geotools

texts:
I need to color polygons in different colors, ensuring that polygons sharing boundaries do not have the same color (the well-known four color theorem).

Is there a way to do that using GeoTools?


texts: Help getting my flex service publicly accessible

texts:
My organization is trying to get our web maps up and running. I have created a simple web map service on our in house ArcGIS server, using the Arcgis flexviewer that is visible locally, works great, and is visible on arcgis.com.

We have a 2 server set up where Server A is our internal server that everyone in the company stores their data on. this is the machine we have our install of ArcGIS 10.1 for server on. It is not public, and for internal staff use only.Server B on the other hand is off site and set up by a company that hosts our website to specifically host our web maps. It is in a DMZ and set up to get through our company firewall. It is just a simple web server running IIS 7 and has the ArcGIS Web Adaptor installed .They, like me, have no experience in setting something like this up so I'm learning as much as I can along the way..

We have our flex map up and running [here] which is great. What is not so great is the error that is popping up when the flex map is opened. It looks like my data, which consists of a a few poly lines, some points and some labels, are not making it over to Server B from Server A. I guess my real question is what do I need to do to get the data from Server A show up in my flex app on server B.

If anyone can shed some light on this that would be great.


texts:
So I have two layers. First one is line layer (columns: id, length, address_id) and atlas is connected to it so when I change objects, map is zooming to next line. Secound is point layer (columns: address_id, city, street, house_number). Both of these layers are related by address_id value and also line layer is intersecting wtih point layer.
I want to display below map, tables in which there will be placed attribbutes from point layer(city, street and house_number) and every time I change atlas object new attributed will be shown in these tables. I was trying to make expressing using attribute function but no success.
Is there any easy way to make it happen. I could insert all information from point to line but i would like to avoid this solution.


texts: Accessing Geoserver plugins through the REST API

texts:
I am writing Python scripts to automate the creation of Geoserver layers. For this the gsconfig module has been invaluable, but I often need to execute operations that are handled by plugins, most often the ImagePyramid one.

I suppose I could use the REST API but I'm not clear how to access the functionalities of installed plugins - in my case ImagePyramid.

Editing to add a more focused reformulation of my question: how do I call on the functionalities of the ImagePyramid plugin using the REST API in Python? More generally, how do I figure it out for any given plugin?