texts:
I am using Geoserver 2.4.4 with a PostGIS database.

I have 2 layers configured, one of them represents some polygons.
I am trying to put a label inside the polygon with the number of the features of the second layer inside the bounding box of the first layer.

Eg: number of cities in a State.

I am wondering if it is possible to achieve this only with SLD styling of the layer.

I am looking at ways to apply Point Stacker function to achieve this, but I can't find a good example or documentation.


texts: Esri ArcSDE 10.3.1 Query Inspire Metadata from Oracle (11g) Enterprise GeoDB via SQL

texts:
As the title says - I want to create a query in Oracle SQL which will read all dates from the geoprocessing history which can be found using ArcCatalog to show the "Description"-tab from a feature. (You need to change the settings of ArcCatalog to use the following "Metadata Style")

If possible I just want to use SQL only - no python or arcpy.

The first approach I used was to query the "cdate" field within "sde.layers" table with the following query:

SELECT Upper(sde.layers.owner
         || '.'
         || sde.layers.table_name)                                    AS
   feature_name,
   ( To_date('1970-01-01', 'YYYY-MM-DD') + sde.layers.cdate / 86400 ) AS
   last_update
FROM   sde.layers
WHERE  sde.layers.owner LIKE 'SDE_V_%'
ORDER  BY last_update DESC;

The query above will list the creation date of every layer within the ArcSDE which is not what I need.

I tried querying the SDE_XML_DOC2 table which should contain the information I described above but wasn't successful because oracle developer always outputs the following error:

Query:

SELECT Xmltype(xml_doc, 871)
FROM   sde.sde_xml_doc2
WHERE  sde_xml_id = 8641;

Error:

LPX-00217 invalid character error

I tried different character sets (1 US7ASCII, 170 WIN1250, 171 WIN1251, 178 WIN1252, 871 UTF-8) within the XMLType function but all of them returned the same error.

I'am stuck at this point and need some help. The general idea is to get the last modified date of all layers within the sde (vector and raster).


texts:
I'm doing this sql query using geoservers sql views:

SELECT kmeans, count(*), ST_BUFFER(ST_CENTROID(ST_Collect(geom)), qrt(count(*))*1000  ) AS geom
FROM (
   SELECT kmeans(ARRAY[ST_X(ST_CENTROID(geom)), ST_Y(ST_CENTROID(geom))], 10) OVER (), geom
      FROM data WHERE foo = %FOO%
)
AS ksub
GROUP BY kmeans
ORDER BY kmeans

This always returns the same regardless which bounding box is requestet via the wms. But I want to use the bounding box in my WHERE clause. I know that the user gets a different result then if there are  more tiles used to dislpay a map but that will be ok for me.

Regards

Daniel


texts:
I need to extract geotiff from SMAP DATA (same without registration)

I found example of gdal_translate command that do this task:

gdal_translate -a_ullr -17334194 7356860 17334194 -7356860 -a_srs epsg:3975 -a_nodata -9999 HDF5:"C:\files\SMAP_L3_SM_P_E_20170711_R15060_001.h5"://Soil_Moisture_Retrieval_Data_AM/soil_moisture sm0000.tif

It's produce geotiff, but it's projection is not correct. There is some gap:

Here is some info that I found:
EPSG:3975 is the Original EASE Grid, while SMAP data is projected onto EASE Grid2.0 (i.e. EPSG:6933)

+proj=cea +lat 0=0 +lon 0=0 +lat 1=30 +x 0=0 +y 0=0 +ellps=WGS84 +datum=WGS84 +units=m"

What is gdal_translate string for setting correct data-projection.

I tried to set EPSG:6933 in QGIS after reprojection to EPSG:3975. Images is shifting to very-very little to more correct position. Something like half-pixel.


texts: Split disconnected polygons into separate features after using clip tool?

texts:
I'm working with a series of forest polygons (geodatabase feature classes), which I'm trying to clip to an area with a specific land use history (marked as "undisturbed" in the screenshot below)in ArcGIS 10.3.1. I'm then hoping to filter out all polygons less than 1 hectare for additional analysis.

My problem is that some of of the forest polygons ("cove_ecogroup" below) overlap the larger land use polygon in multiple areas, which leaves me with some disconnected features in my output. This is fine, except that it makes my 1-hectare filter impossible because these shapes are being stored as larger than they actually are because they retain the FID of the original feature instead of splitting into multiple new features (See image 2).

Is there a way to split these into multiple features so that I can filter out everything under 1 hectare? This analysis covers hundreds of thousands of acres, and so I won't be able to edit each new polygon individually.


texts: QGIS: QgsMapCanvas rotate vs setRotation

texts:
Question regarding QgsMapCanvas using python for QGIS :

what is the difference between  the following two methods?

myQgsMapCanvas.rotate(10.0)
  myQgsMapCanvas.setRotation(10.0)

(myQgsMapCanvas is of type QgsMapCanvas)

My understanding after trial and error:
When using myQgsMapCanvas.rotate(10.0), the map is rotated relative 10.0 deg and the rotation widget at the the QGIS status bar does not change

When using myQgsMapCanvas.setRotation(10.0)  the map is rotated to absolute 10.0 deg and the rotation widget at the the QGIS status bar changes to 10.0 deg

Is this correctly understood ?

I use QGIS 3.2.1-Bonn and QT 5.9.2


texts: Relocating & resizing complex shapefile with OpenSource Software

texts:
(please excuse my poor english)

I do have a polygone shapefile of an archaelogical site with elevation data of 10cm steps and I dont get hands on the raster data as it was collected by the drone.

My problem is, that I need to move and resize the polygones, due to even though the shape say's it is a DHDN/GK3, the lower-left corner of the collected data is 0,0 and the x & y extends seam to be decimenter relative to this origin.

Is there any useful way to solve the problem with QGIS + tools (using Ubuntu Linux) ?

My so far best thought is....

a) CSV-Geometry export via MMQGIS (temp-nodes.csv, temp-attributes.csv)
   The CSV counts ~8.500.000 nodes
b) Suggest origin from so far data (kept as variable)
c) pipe node-csv through script to convert decimeter to meter,
   to get the new delta x/y from the origin
d) add the new delta x/y to the origin and save
e) CSV-Geometry improt via MMQGIS

Is there any way to solve this easier?
It's been years since I have written scripts.

Chris


texts:
Last week we used Collector for ArcGIS (version 10.3 on an iPad Air, AT&T) with our ArcGIS Online organization account for field data collection. We experienced loss of data, sometimes everything including the point feature itself, other times photo attachments. Since the collection was at or near urban areas and cellular signal strength was good in general, we didn't use offline mode and it was all done real time.

On day 2 of the collection there were a few instances with multi-photo points when the “submit” button was tapped it would take the App a while to process then a submission failed message would appear, offering two options – “Cancel” and “Retry”. Initially “Retry” was selected but the same failure message would appear again and again, so eventually the field crew clicked on “Cancel”. At that point the collected feature as well as attributes and photos were visible on the App, so they moved on to the next site.

When they were back at the hotel (on Wi-Fi) they realized 2 points they collected during the day were completely gone, and 3 other points were missing all the photo attachments. They also noticed 2 points from day 1 were completely gone (they couldn't recall with certainty whether the “failed to submit” message was encountered with them). Not knowing how to avoid more data loss, they went with paper and pen for the remaining of the collection and re-collection of the lost points.

Once they were back in the office we did more testing. While the cellular signal strength was at least 3 bars out of 5, it was pretty easy to duplicate the issue by taking more than 3 or 4 photos and submit.

Here are our questions:

Is our instinct correct that our issue could be caused by the data volume submitted each time?
We like the convenience and ease of set up with this solution. However based on our first project the reliability is too low for serious application. I can’t imagine we are the only one needing multiple photos for each collected features. Has data loss been a common issue?
What is everyone else doing to make this solution work? I assume the reliability will improve if each submission is limited to one photo, but it’s such an interruption to the workflow.
When collecting in “online mode”, is there a way to explicitly back up the data collected onto the mobile device?
We are thinking maybe offline mode can be more reliable and will do some tests. But from experience, will syncing (Wi-Fi or cellular) be problematic if a collected feature contains a large number of photo attachments?
Any recommendation on more reliable solutions that works with ArcGIS Online?
How about complete on device Apps?

Update: after talking with an ESRI technical support back and forth, he started being able to duplicate the problem on his end, and has registered this issue as a bug. We did perform more tests of our own. It seems that one photo per submission/update is problem free, so while on cellular only if multi photos are needed for a collection point you'll have to update after each photo is taken. Offline mode seems to be working fine.


texts: Seeking options for Spatial ETL (Extract, Transform, Load)?

texts:
I am interested in the pros and cons of various spatial ETL (extract, transform, load) tools.  If you have used the items listed here (or add your own), I seek your opinions and experiences.  In particular I would like to see usability comparisons of:

Safe Software FME (proprietary)
ArcGIS Data Interoperability Extension (proprietary, subset of FME features)
Talend Spatial extension (previously Spatial Data Integrator) (open source)
GeoKettle (open source)
ogr2ogr (I don't know that it is truly ETL, but more of a file translation toolset)
Any additional lot listed

There is no need to give a review of ALL software mentioned.  If you are experienced with even one then that will be very beneficial in making a decision about which direction to go.

Example: I am looking to create a schema conversion function that will allow me to select the input layer, create a translation, and output to a new, pre-defined schema.  Optimally, after creating the translation script, I would like to have an interactive form where I can "map" fields in my input layer to the output layer (ie- The output layer will have a field called "Address", what is it called in the input layer?)

Some were mentioned in the Q&A at What tools are available for uploading gis data to a database?

And here are a couple of related articles that I found.

http://fuzzytolerance.info/code/open-source-etl-tools/
http://geobiwiki.spagoworld.org/xwiki/bin/view/Main/Talend+GeoKettle
http://slashgeo.org/2011/07/25/Open-Source-ETL-Tools-GeoKettle-20-and-Spatial-Data-Integrator-420-Released


texts: Converting floating point raster to vector

texts:
I'm trying to convert a raster layer to vector but I'm really struggling.
First it says that I have floating points rather than integers so I have I have seen you can use raster calculator to multiply the numbers to make them whole. Multiplying the points did not work, it would not let me do it - I do not know the error code right now for it but I can have a look. What I am trying to do is extract the raster from the vector but have ranges for the classifications.
I am using ArcGIS Pro but if your help requires me to use QGIS I will be happy to try. I also have seen that some kinda solutions require me to code but I have no experience with coding.
Well I am trying to convert to a vector as I want to understand the percentage of surface area which is an ideal climatic location for a species. I do not want a polygon around each pixel and want to add ranges or categories and have each pixel fit into a range so then I can have larger polygons. Or if that is not possible to extract from the raster a certain range so I can see the different areas.


texts: Programmatically Harvest and Disapprove Metadata - GeoPortal Server

texts:
Is it possible to synchronize an existing metadata catalog in esri's GeoPortal Server? I have access to the server, and have a script to add an XML file to our web accessible folder. After this script runs, I manually log in to the Geoportal site and synchronize content.

Is there any way to run this synchronization without using the web interface? Ideally it would be run in the command line on server, but I could look into (somehow) calling this from Java as well, as may be described in the GeoPortal documentation. Do I need to get into extending the web harvester or is there an easier way to do this?

Finally, is there a way to programmatically set the status of an existing record by GUID (i.e. approved/disapproved)?


texts: How much a decimal degree is equal to in latitude 40N?

texts: I would like to learn how much a decimal degree equals to 40N latitude as units of distance. I have this table (http://en.wikipedia.org/wiki/Decimal_degrees) in mind, but I am interested for the latitude 40N.


texts:
I'm working with openlayers3 and I wanted to load a GeoJSON file from my local ../assets/data/countries.geojson. Nothing is wrong but it didn't appear on the screen.

This is my code:

<!doctype html>

    <title>Vector Examples</title>
    <link rel="stylesheet" href="http://openlayers.org/en/v3.13.0/css/ol.css" type="text/css">


texts: How to hide some parts of an image layer with OpenLayers?

texts:
My Problem

I'm currently displaying a layer image above a Google Maps base layer. I would like, when the user check/uncheck some checkboxes, that some parts of this image layer disappear, and the base layer appears instead these parts.
But I want the other parts still visible with the same opacity.

Explanations

The image parts to hide are polygons. Those polygons are already generated as layers and are above my image layer. I could put a black background color to hide the wanted image parts, but it wouldn't display the base layer instead...
I'm stuck.

Solutions I thought about

I would like you to say me if that's possible and, if you can, how to do it.

I'm not asking for some code, I can do it myself, I'm just kind of a newbie at OpenLayers and I'd like someone to help me to understand how it works. Plus, I have to add functionnalities to an already written code so it's not very easy.

Here are data I have, I removed the ones which are not useful (this is only an example, empty data are usually filled with correct values).

The images :

images[id] = {
    name: "Image",
    image_date: '2012/08/25',
    image_path : 'static/path/to/image',
    width : '250',
    height : '250',
    georef_north : '',
    georef_south : '',
    georef_east : '',
    georef_west : ''
};

And the polygons :

polygons[id] = {
    name:'',
    longitude:'',
    latitude:''
};

I'm getting them like that via WFS protocol like this :

layers['blockLayer'][i] = new OpenLayers.Layer.Vector("Block", {
        strategies: [new OpenLayers.Strategy.BBOX()],
        filter: VMMS.GIS.GetBlockFilter(i),
        protocol: new OpenLayers.Protocol.WFS({
            url:  "../vine/jsGeoserver",
            version: "1.1.0",
            featureType: "block",
            featureNS : urlGeoserver + dataStore ,
            featurePrefix : dataStore ,
            srsName: "EPSG:900913",
            //srsName: "EPSG:4326",
            geometryName: "block_geom",
            schema : urlGeoserver +"wfs/DescribeFeatureType?version=1.1.0&typename=" + dataStore +":block"

        }),
        projection : proj_spherMerc,
        //projection : proj_WGS84,
        styleMap : StyleBlockWithOpacity
    });

First Solution

  Why wouldn't we put the image layer in each polygon layer, and display only the layers we want to display ?

Would that work ? I'm not actually sure about that and I don't want to get into it if it's not realizable. Because I don't know if I can display an image IN a layer and not entire, and it has to be well placed in the map and in the layer... Though it seems to me it's doable.

Second Solution

  Why wouldn't we display the base layer again above the image to hide some parts ? The image would still be here, but as the base layer would be on top of it, we wouldn't see it.

This solution seems to me even less doable, or, at least, very heavy.

An example image, it's what I'm currently displaying : There's the google maps base layer, then the image blue and green, and all forms are polygons layers above the image.
Screenshot

Thank you all for your answers !


texts:
We would like to render a polygon on our MapServer.  Everything works wonderfully on our localhost.

Now we would like to send the application via external server. Everything works fine on the external server as well. Except that the polygon appears on the map. The polygon is a TabFile. All paths have been updated, adapted and checked several times.

In the network analysis comes with the "map".exe this error: "No query information to decode. QUERY_STRING is set, but empty."

Someone an idea what else would have to be adjusted?


texts:
I would like to have a pop-up menu when right-clicking the MapCanvas
(the identifyFeatures tool does it...)

Could anybody give me a starting point on how to achieve this?


texts: Force tiles to load higher resolution in Openlayers

texts:
I've got some offline pbf files that I load into Openlayers through a HttpServer. These are static vector tiles of a certain area. Now I'm working on making sure that these are visible at the zoomlevel it needs to be.

From the data, the maximum zoom level is 14. I'm trying to view it at 18. Which has worked by loading the map at zoom level 14 before switching to 18 (with a very hacky setTimeout). After moving too far on the map towards the next tile, it won't show it anymore because it requests tiles at zoomLevel 18. Obviously this isn't good practice.

How else would one do this? The basic behaviour is that at Zoom 18 I want to display tiles from Zoom 14.

code:

  var baseLayer =  new ol.layer.VectorTile({
            declutter: true,
            renderMode: 'vector',
            preload: Infinity,
            source: new ol.source.VectorTile({
                format: new ol.format.MVT(),
                url: this.tilesUrl,
                overlaps: false,
                tileGrid: new ol.tilegrid.TileGrid({
                    extent: ol.proj.transformExtent(this.route.extent, 'EPSG:4326', 'EPSG:3857')
                })
            }),
            style: this.style
        });

What I've looked into so far:

Set 'Preload' on the source to true.
Setting a TileGrid (I couldn't get this one to work properly)
Setting custom resolutions myself and forcing the resolution to zoom 14 and below (this just simply didn't work, it won't load any data at all trying this)


texts: Seeking tools for uploading GIS data to database?

texts:
I would like a nice graphical tool that allows me to take standard gis data such as shp files and kml files and upload them to a database such as PostGIS, MySql or Oracle.

Are there any such tools?


texts: Bulk Exporting Shapefiles Into Microsoft SQL Server using QGIS or PostGIS

texts:
I have found two popularly suggested ways to preform this task:

ogr2ogr - just take a shapefile and use ogr2ogr to put a shapefile in to a sql server database.
This would work for a few shapefiles but is not a solution to importing LARGE amounts of shapefiles of many varieties.

Using QGIS, export into SQlite - I have had this one proposed to me a few times. Really this only adds the steps necessary to preform this to many many shapefiles.

And this third task I think may be the winner:

Exporting shapefiles into Postgis, PostgreSQL, then perform a batch database conversion from PostgreSQL into Microsoft SQL Server 2012 - After breaking down and buying software called FullConvert which claims to be able to connect the the PostgreSQL server and then copy tables into MS SQL Server in bulk.

The last method I think is the closest I have gotten to being able to preform this task. I'm going into the office tomorrow to preform this on the server locally because it fails every time I try to perform this task over OpenVPN.


texts: Raster Combine using ArcGIS Spatial Analyst?

texts:
I am trying to combine a whole bunch of rasters with ArcGIS Spatial Analyst.  Seems to be doing what I want EXCEPT the calculated output is only the area where all 4 rasters overlap.

I am new to this, and don't know what I'm doing wrong!

Using 9.3.1


texts: Linux based QGIS File Geodatabase access

texts: I have looked for a "recipe" to get QGIS to read FGDB but cannot find anything. I'm fairly new to Linux and need a step by step procedure or somewhere to find the steps. I'm running QGIS 2.14.3 Essen on Mint.


texts:
This is similar to this rather old post for a much older version of QGIS, but I'm still getting the same behaviour: Defining new custom projections in QGIS?

Basically I'm defining a new Custom Projection in QGIS using this proj4 string from epsg.io (http://epsg.io/26332-1167):

+proj=utm +zone=32 +ellps=clrk80 +towgs84=-81,-84,115,0,0,0,0 +units=m +no_defs

I'm entering it as per the screenshot below. I click OK.

Then I reopen Custom Projections straight after, and all the proj4 params have completely changed:

Why is this happening?

Are the params in the second screenshot equivalent to those specified in the first?

I'm on QGIS version 3.10.1, the screenshot below shows the other gubbins if useful:


texts: Calculating point GPS location in picture

texts:
I have camera location (GPS coordinates) and need to calculate some points (pixels, actually it’s some detected objects points) location (GPS coordinates) that are seen from the camera. I found in the camera manual it’s viewing angle (Opposite Angle 156°, Horizontal 122°, Vertical 89°) but no other parameters that would seem to be useful. Is it possible to calculate the location of some pixels? Could you give me some hints?
I found some examples of calculating distance from the camera, but I think it’s a different problem.


texts:
I have a KMZ file of north-eastern Thailand which marks the boundaries of a larger district and it's sub-districts. Still new to GIS and in the process of getting familiar with QGIS, I wonder can anyone advise me how can I "extract" the boundaries of each of these sub-districts from that KMZ?

E.g. are there any QGIS plugins which could be used for this task?


texts:
I'm newbie in coding and I'm trying to modify some downloaded GitHub code to make customized application and display my data.

AddLayer function is making me desperate.
I'm trying to add data from MapBox account by changing  the part of code with my Mapbox username and project id:
     L.mapbox.tileLayer('myMBid.n7i50n7g'),
      L.mapbox.gridLayer('myMBid.n7i50n7g'),

Please share your coding knowledge with me and teach me how can I modify code to be able to load data from MapBox?

Or how can I, inside this coding, insert option to load geojson or topojson files from application folder?

Full Code:

<pre>
var map = L.map('map',{maxZoom:17,minZoom:1}).setView([45.2701486, 18.5586462], 13);

L.mapbox.accessToken = 'pk.myMBpublickey';

var baseLayer = new L.mapbox.tileLayer('myMBid.hl6099hm').addTo(map);

// ADD LAYER CONTROLLER
var ui = document.getElementById('layerControls');

//ADD FIRST LAYER
addLayer(
  L.mapbox.tileLayer('myMBid.j2rlctfi'),
  L.mapbox.gridLayer('myMBid.j2rlctfi'),
  'Administrativne granice',
  1
);

//ADD SECOND LAYER
addLayer(
  L.mapbox.tileLayer('myMBid.n7i50n7g'),
  L.mapbox.gridLayer('myMBid.n7i50n7g'),
  'Soil Types',
  3
);

function addLayer(layer, gridlayer, name, zIndex) {
  layer.setZIndex(zIndex);
  var gridControl = L.mapbox.gridControl(gridlayer, {position: 'bottomleft'}).addTo(map);
  // Create a simple layer switcher that toggles layers on
  // and off.
  var link = document.createElement('a');
  link.href = '#';
  link.className = 'btn btn-primary btn-sm';
  link.type = 'button';
  link.innerHTML = name;
  link.onclick = function(e) {
    e.preventDefault();
    e.stopPropagation();
    if (map.hasLayer(layer)) {
      map.removeLayer(layer);
      map.removeControl(gridControl);
      this.className = 'btn btn-primary btn-sm';
    } else {
      map.addLayer(layer);
      map.addLayer(gridlayer);
      map.addControl(gridControl);
      this.className = 'active btn btn-primary btn-sm';
    }
  };
  ui.appendChild(link);
};

// ADD THE REFERENCE OVERLAY
var topPane = L.DomUtil.create('div', 'leaflet-top-pane', map.getPanes().mapPane);
var topLayer = new L.mapbox.tileLayer('landplanner.hl60jemk', {
  maxZoom: 18
}).addTo(map);
topPane.appendChild(topLayer.getContainer());
topLayer.setZIndex(7);
</pre>


texts: ST_Multi returning GeometryCollection in PostGIS?

texts:
I am working in Postgres 9.6. I have a table called builtup_areas:

                                          Table "public.builtup_areas"
    Column    |             Type             | Collation | Nullable |                    Default
--------------+------------------------------+-----------+----------+------------------------------------------------
 ogc_fid      | integer                      |           | not null | nextval('builtup_areas_ogc_fid_seq'::regclass)
 objectid     | bigint                       |           |          |
 wkb_geometry | geometry(MultiPolygon,27700) |           |          |

And a table called local_authorities:

                                           Table "public.local_authorities"
     Column     |           Type           | Collation | Nullable |                      Default
----------------+--------------------------+-----------+----------+----------------------------------------------------
 ogc_fid        | integer                  |           | not null | nextval('local_authorities_ogc_fid_seq'::regclass)
 wkb_geometry   | geometry(Geometry,27700) |           |          |
 canonical_name | character varying        |           |          |

I am trying to create a new table that contains the ST_Difference between the local authority's boundary, and the built-up areas, for each local authority. I have made a new table to hold this, with a MultiPolygon column:

CREATE TABLE las_without_buas(gid serial primary key, geom geometry(MultiPolygon,27700), county varchar(1000));

This is my code:

INSERT INTO las_without_buas(county, geom)
  SELECT
    county,
    ST_Multi(ST_Difference(la.geom, bua.wkb_geometry)) AS geom
  FROM
    (SELECT canonical_name AS county, wkb_geometry AS geom FROM local_authorities) la,
  builtup_areas bua;

But when I run it, I get this error:

ERROR:  Geometry type (GeometryCollection) does not match column type (MultiPolygon)

I don't understand why it's trying to insert a GeometryCollection when I've specifically used ST_Multi to return a MultiPolygon. Any suggestions?


texts: Can't clear markers from MarkerClusterGroup

texts:
i have problems removing all markers befor adding new ones.

  var L_Emergency_Access_Points = L.markerClusterGroup();
  var EapMarkers = L.geoJson(null,merge_options(globalOverlayOptions,
               { name:          "EapMarkers"
                ,onEachFeature: function (feature, layer) {
                                   layer.bindPopup("hallo")
                                }
               })
            );

... snip ...

map.on("moveend") {
   bbox = map.getBounds().getSouth().toFixed(5) +","+  map.getBounds().getWest().toFixed(5) +","
        + map.getBounds().getNorth().toFixed(5) +","+  map.getBounds().getEast().toFixed(5);

    $.ajax({
           url: geosWFS,
           async: false,
           data: {
                   service : "WFS",
                   version : "1.1.0",
                   request : "GetFeature",
              outputFormat : "application/json",
                   SrsName : "EPSG:4326",
                   typeName: "osm:Emergency Access Points",
                       bbox: bbox
           },
           datatype: "json",
            success: function(data,text,xhr) {
                       console.log(myLocalName, "success");

                       EapMarkers.clearLayers(); // problem ????
                       map.removeLayer(L_Emergency_Access_Points);

                       $(data.features).each(function(key, data) {
                          EapMarkers.addData(data);
                       });

                       L_Emergency_Access_Points.addLayer(EapMarkers);
                       map.addLayer(L_Emergency_Access_Points);
                    }
        });
});

after loading page:

after litte panning to trigger "on moveend":

reload data is triggered by "on moveend" and it is working fine. But it does not clear the old data. every pan the new markers are added without removing the old ones.

btw: can't add more images :( your restrictions to new users are not helpfull. :(

Any idea?

EDIT: Got it by replacing success block:

                     success: function(data,text,xhr) {
                             console.log(myLocalName, "success");

                             L_Emergency_Access_Points.clearLayers();
                             map.removeLayer(L_Emergency_Access_Points);

                             var markerList = [];
                             $(data.features).each(function(key, data) {
                                var marker = L.marker(L.latLng(data.geometry.coordinates[1],
                                                               data.geometry.coordinates[0]),
                                            { title: "hallo" } );
                                marker.bindPopup("hallo");
                                markerList.push(marker);
                             });

                             L_Emergency_Access_Points.addLayers(markerList);
                             map.addLayer(L_Emergency_Access_Points);
                          }


texts: arcpy to handle topology error inspector

texts:
I am now want to do topology check use python, but I don't know whether arcpy can realize the function in error inspector. May I call error inspector and get the error records(something like attribute table) using python(when arcmap is open)?
Thanks!


texts: 3D map view in QGIS with .grd format bathemetry data

texts:
I am attempting to use the 3D map feature available in QGIS3 for viewing bathymetry data. These data are in .grd format, but I’ve also tried using data in GeoTiFF format with no luck.

The 3D map comes up blank, and no amount of scrolling brings the panel into view. I’ve checked that this feature does, in fact, work with .hgt (SRTM) elevation data.

I’m not super familiar with the differences between different data formats (.tif, .grd, .hgt, etc...) but my thinking is that these data are either missing some information or need to be converted to .hgt.

I do know that conversion to .hgt is possible, but it has some drawbacks that I’d like to avoid.

Has anyone encountered this problem and know of a possible solution?

EDIT: Here are some data


texts: combine rasters into single raster with multiple field values in R

texts:
I have several raster brick objects that I would like to combine into one raster with multiple field values for each cell, but I am not sure how to do this in R.

For example, let's say I have two layers in a 2*2 raster brick.

r <-raster(nrow=2,ncol=2)
r[] <- round(runif(ncell(r))* 10,0)
s <-raster(nrow=2,ncol=2)
s[] <- round(runif(ncell(r))* 10,0)
t = brick(r,s)

I would like to combine the values for r and s into a single raster with multiple values (at the object@data@attributes slot, I think) associated with each cell.

So, in the above example, if cell 1 = 1 in r and cell 1 = 2 in s, then the output that I desire is that the resulting raster contain two fields at @data@attributes, the first field = 1 and the second field = 2.


texts: Creating a line from lots of easting/northing co-ordinates

texts: I have hundreds of easting northing points which I need to join as a single line to show a route.  How do I import them and get QGIS to automatically turn them into a line?


texts: Processing runalg gdalogr Error: Algorithm not found

texts:
I can not use the "gdalogr"-Processing in  QGIS 3.0 (2.99)
QGIS code revision 1182816

import processing
processing.runalg('gdalogr:convertformat','', 0, '' ,'')
Error: Algorithm not found

Does anyone know whats going  wrong?


texts: Launching one widget from a button inside another widget

texts:
Suppose I have created a button inside a custom InfoTemplate (aka 'popup') widget in my Web AppBuilder app:

 var editButton = new Button({
     label: "Edit Schedule",
     onClick: function () {
         editButton.setDisabled(true); //disable the edit button on first click
     }
 });

 editButton.startup();

What I want is to launch another widget called 'TaskManager' when the onClick event is triggered.
It seems like I need to use a library class called WidgetManager but I can find any practical example for my scenario. It seems like I need to use something like this:

 openWidget(widget)

The openWidget method requires the Widget ID, which I am having trouble locating in my app (I know the name is TaskManager from the manifest.json file but can't find the id attribute).
I've also seen example where a function is created, I assume inside the InfoTemplate widget, that will point to the TaskManager widget and can be called inside the click event:

   _openTaskManager: function () {
       var taskManager, sbc;
       var widgetCfg = this._getWidgetConfig('CustomInfoTemplate');
       if (widgetCfg) {
           sbc = WidgetManager.getInstance().getWidgetsByName('TaskManager')[0];
           sbc._resizeToMax();
           sbc.setOpenedIds([widgetCfg.id]);
       }
  },

then => _openTaskManager();

Am I on the right track? Is there an easy way to do this or am I way off? Any suggestions?


texts: Cannot add multiple basemaps at once using Python console of QGIS

texts:
I am using Qgis 3.8. I want to download a bunch of basemaps to have easy access to in the future. I tried using python console with URl: https://raw.githubusercontent.com/klakar/QGIS_resources/master/collections/Geosupportsystem/python/qgis_basemaps.py.
I get this error message:

  SyntaxError: invalid syntax.

What am I doing wrong?


texts: Relationship classes and feature classes with the same name

texts:
I've noticed that it is possible (at least in ArcGIS 10.1) having within a Personal geodatabase, a relationship class and a feature class/table with the same name!. This is my example:

The above example was created by using ArcCatalog manually, but now I want to do some scripting in Python.

The problem I'm having is that after using arcpy.Describe(path_to_dataset).children, in which I can easily list the existing relationship classes by using datatype == "RelationshipClass". Then, I'm unable to access the rc individually on this way:

arcpy.Describe(path_to_rs) because ArcGIS picks up by default the fc always, not the rc, and both have the same path.

Any ideas?


texts: How to fix this "Sorry  is not a valid option" in GRASS-Ubuntu using r.water.outlet?

texts:
I am working with this .sh code (a while loop, extracted from: http://www.surfaces.co.il/) to create watersheds from drainage raster and outlet point of every watershed, but I obtain this error: "Sorry <1> is not a valid option", "Sorry <2> is not a valid option" and so on.

i=0   #an iterator to give consecutive names to the new watersheds
while read X Y; do
    i=$(( i + 1 ))
    r.water.outlet drainage=fdir basin=wshed$i easting=$X northing=$Y --overwrite
    echo "$X"
    echo "$Y"
    #r.null -q wshed_$i setnull=0
    r.to.vect -s in=wshed$i out=wshed$i feature=area --overwrite
    #g.remove rast=wshed_$i
done < cross_points.txt
echo "Created $i catchments"


texts: QgsMapCanvas.setExtent() has no effect when loading and rendering many layers

texts:
The following script is intended to load a pile of vector layers (50 times the same) and after loading to set map canvas center and scale:

from qgis.core import QgsProject, QgsVectorLayer, QgsDataSourceUri, QgsPointXY
from qgis.utils import iface

canvas = iface.mapCanvas()
center = QgsPointXY(4428590, 5495530)
root = QgsProject.instance().layerTreeRoot()
current_group = root.addGroup('Test Group')

layers = 50 * ['routing_network']

u = QgsDataSourceUri()
u.setSrid('31468')
u.setDatabase('c:/tmp/routing_slow.sqlite')
u.setKeyColumn('ogc_fid')
for lyr in layers:
    u.setDataSource('', lyr, 'geometry')
    l = QgsVectorLayer(u.uri(), 'test layer', 'spatialite')
    QgsProject.instance().addMapLayer(l, False)
    current_group.addLayer(l)

canvas.setCenter(center)
canvas.zoomScale(2000)
print('new extent set')

When the for loop has run through, the scripts prints out new extent set immediatly, but leaves map canvas extent to the full extent of the layers:

This is how it should look, when setting things on the console after rendering is completed:

Since rendering the layers takes some (~10) seconds, can it be that the canvas extent is set to the layers extent after rendering overriding the setCenter() and zoomScale() at the very end of the script?

Weird observation: when doing a iface.mapCanvas().zoomScale(50000) (scale in full extent is 1:192000) during rendering rendering continues at scale 1:50.000.

How to deal with this? Is there a way of setting map canvas extent only when rendering is complete?


texts:
I created some image overviews for a large (1.2 Gb) topo relief image.

gdaladdo -r cubic -ro --config COMPRESS_OVERVIEW JPEG --config PHOTOMETRIC_OVERVIEW YCBCR --config BIGTIFF_OVERVIEW YES relevo_srtm_br_30m.tif 2 4 8 16 32

But as a result I'm getting some stripes on my overviews

However, when I zoom in enough to look at the original data, the stripes disappear.

Is this beeing caused by the cubic interpolation?


texts: Problem with WFS-T save strategy in Openlayers

texts:
In my web map application, I have added the editable wfs vector layer for feature editing. But in save strategy method, Insert and Update Geometry is working fine.But Update Attributes is not saved in transaction. Did I miss anything to include while adding editable layer? Please suggest me to fix this problem.

var saveStrategy;
function init(){
format = 'image/png';

bounds = new OpenLayers.Bounds(
            -288161.9289062237, -435.65815425360597,
            949406.755566389, 1300000
        );
var options = {
controls: [],
maxExtent: bounds,
maxResolution: 5079.826789665053,
projection: "EPSG:404000",
units: 'm'
};

 osMap = new OpenLayers.Map('map', options);

//Setup Base layer
tiled = new OpenLayers.Layer.WMS(
            "Oxford_Map", "http://xxx.xx.x.xx:8080/geoserver/cite/wms",
            {
                LAYERS: 'Oxford_Map',
                STYLES: '',
                format: format,
                tiled: true,
                tilesOrigin: osMap.maxExtent.left + ',' + osMap.maxExtent.bottom
            },
            {
                buffer: 0,
                displayOutsideMaxExtent: true,
                isBaseLayer: true,
                yx: { 'EPSG:404000': false }
            }
        );

untiled = new OpenLayers.Layer.WMS(
            "Oxford_Map", "http://xxx.xx.x.xx:8080/geoserver/cite/wms",
            {
                LAYERS: 'Oxford_Map',
                STYLES: '',
                format: format
            },
            {
                singleTile: true,
                ratio: 1,
                isBaseLayer: true,
                yx: { 'EPSG:404000': false }
            }
        );

 osMap.addLayers([untiled, tiled]);

 //set up a save strategy
saveStrategy = new OpenLayers.Strategy.Save();
saveStrategy.events.register("success", '', ChangesSuccess);
saveStrategy.events.register("failure", '', ChangesFailed);

wfs = new OpenLayers.Layer.Vector("ODC_Private_Graphics_Data", {
    strategies: [new OpenLayers.Strategy.BBOX(), saveStrategy],
    projection: 'EPSG:27700',
    protocol: new OpenLayers.Protocol.WFS({
        version: '1.1.0',
        url: 'http://xxx.xxx.x.xx:8080/geoserver/wfs',
        featurePrefix: 'cite',
        featureNS: 'http://www.opengeospatial.net/cite',
        featureType: 'ODC_Private_Graphics_Data',
        geometryName: 'the_geom',
        extractAttributes: true,
        srsName: 'EPSG:27700',
        isBaseLayer: false,
        visibility: true,
        schema: 'http://xxx.xxx.x.xx:8080/geoserver/wfs/DescribeFeatureType?version=1.1.0&;typename=cite:ODC_Private_Graphics_Data',
    });

osMap.addLayers([wfs]);
}

    //This function will be invoked when I click on save image in web map application
 function saveFeature() {
  for (var ss = 0; ss < wfs.selectedFeatures.length; ss++) {
      var selFeat = wfs.selectedFeatures[ss];

      selFeat.attributes.FeatureName = 'Test';
       selFeat.state = OpenLayers.State.UPDATE;
  }

  saveStrategy.save();
 }

function delFeature() {
if (wfs.selectedFeatures.length > 0) {
for (var dd = 0; dd < wfs.selectedFeatures.length; dd++) {
     var selFeat = wfs.selectedFeatures[dd];
     if (selFeat.fid == undefined) {
         this.wfs.destroyFeatures([selFeat]);
        } else {
        selFeat.state = OpenLayers.State.DELETE;
        this.wfs.events.triggerEvent("afterfeaturemodified", { feature: selFeat });
        selFeat.renderIntent = "select";
        this.wfs.drawFeature(selFeat);
        }
    }
}

 function ChangesSuccess() {
  alert('Done');
 }

 function ChangesFailed() {
  alert('Failed');
 }


texts: Fixed buffers and interior buffers are not working for a polygon that contains a hole

texts:
I'm following the online qgis manual available on qgis website. I'm doing the step 21.11.3. (http://docs.qgis.org/2.14/en/docs/training_manual/answers/answers.html#complete-analysis-3).

To make it simple : I have a raster converted to vector, so now it is a polygon (several polygons). I'm doing an interior buffer exactly how it is described in the manual. I didn't obtain the same result. I can see the interior buffer just in some polygons, not in all of them ->

The interior buffers appear just for those polygons, that don't contain a hole. I use qgis 2.18.3.

This is how i created the interior buffer ->

In green you see the input layer : suitable_terrain, in blue the output layer, the interior buffer.

It is not a problem of a projection. I use WGS84 UTM Zone 34S where the units are in meters. I did also the Check validity and i obtain the same result.


texts: What techniques exist for clustering road routes?

texts:
The typical clustering (by proximity) techniques that I've seen - hotspots or marker counts - relate to single point data. Different challenges exist for deciding the significance of proximity in multi point line_strings e.g. road routes:

many routes may start together but finish separately
only a single intersection is shared vs. parallel routes

What techniques exist for grouping a number of lines?


texts: Using a geo-referenced image as an overlay to display markers on floorplan

texts:
I've created a tool based on the Google overlay tiler. This allows our users to upload a floor plan, adjust the position etc. using Google maps. Nothing fancy.

We have a server-side job which creates the geo-referenced image using gdal_translate. The job then slices it up using gdal2tiles.py and then saves to s3.

This all works nicely. The idea is that they can use this to correctly position their wireless access points etc. The user can zoom in and out. Again, all is well. Something like this:

The problem I have is that they actually need to see the floor plan in it's original format - in full screen. Not rotated as the image shows.

I don't really need the map at this point - I just the image with the markers in an overlay as so:

I obviously need to rotate the map so the floor-plan is aligned correctly horizontally.

How I can use the data in the geo-referenced image etc. to calculate the scale for the map and the rotation I need?


texts: Is Transformation Required for empty feature classes?

texts:
I am trying to change the projection system of 100+ feature classes from projected to WGS84 using Python.

I am exporting the geodatabase to xml (Schema only)
and changing the projection parameters in xml and then creating a geodatabase from that new xml.
Normally, in ArcGIS Desktop, if I have the data and I change the projection then I have to apply projection transformation as well.

If I do not have any data inside a geodatabase then do I need to apply transformation?


texts: Creating If Then statement in Field Calculator between 2 values?

texts:
I have a FC with a column for NEAR_ANGLE from the output of the NEAR tool. The angle is displayed in degrees (i.e. -154.67687) which is fine. I am trying to create a column, HEADING, that reads this degree value and converts it to either a "North, North-east, South, South-west... etc. whichever direction the near object happens to be from my point.

i.e. If value is between 0 and 90 then East, or if value is between -180 and -90 then South-west

I've searched and tried several examples but none work for me.

Ideas?


texts: How to interpret the OSGeo4w software version under Advanced Install option?

texts:
When choosing to install QGIS Desktop using the OSGeo4w advanced install option, I noticed that there are what looks like ambiguous software versions to choose.

For example, currently QGIS Desktop is:
1.8.0-4
or if you click on the version, it changes to
1.8.0-3

and under the libraries to choose:
OGR FileGDB Driver (Is this the ESRI API)
1.9.2-1
or if you click on the version, it changes to
1.9.0-1

What does the number after the hyphen refer to?  It seems like it could install any number of versions, but likely it has a specific meaning.  Thanks!

(I'm going to install QGIS Desktop with FGDB support, and I think GDAL 1.9.0 is the one that's been proven to work with QGIS)


texts:
Can anyone provide an example of how one would use the (pre-defined) function

Function HttpPost( ByVal strURL As String
                 , ByVal strHeaders As String
                 , ByVal strOptional As String
                 , ByVal bAuthenticate As Logical
                 , strContentType As String
                 , ErrorMsg As String) As CString

in MapBasic?

My POST message needs to send 2 parameters, e.g. findthis and bbox, to a php script sitting on our server which will use the parameters in an SQL query and then returns json-formatted results.

All I've been able to find is sample code within my MapBasic installation folder, however this only uses the HttpGET function and no examples of HttpPOST.

I have not found documentation.

I have posted on MapInfo-L as well where I was able to get a response with HttpGET but then not any further.

I have not yet emailed Pitney Bowes support.

Note: need to do this in 32bit MapInfo (i.e. up to v15), the application I need this for hasn't been converted to 64bit yet


texts: Is there a way to sort the QGIS attribute table by date modified?

texts: I was wondering if there was a way to sort the attribute table in QGIS Pi by date modified. I have been digitizing a large area and I need to delete the last 50-100 features that were digitized. My current table has roughly 3,000 features on it and I really don't want to search through the map or attribute table to find them one by one.


texts: Getting google static map at inbetween zoom levels that google maps site does?

texts:
I'm not sure if this is new in the Google maps site but you can use the scroll bar to zoom in-between the standard zoom levels.  The plus and minus buttons still jump a full zoom level and there doesn't appear to be any change in the raster in the in-between zoom levels although some road labels change (are they vector then?)

Can I get a static map at the in between zoom level as I often want it a bit more zoomed in but not all the way to the next level?

Also is the in-between zoom available in the JavaScript API?

It doesn't look like this is possible from the Google documentation but that's not always very up to date.


texts: Can anybody help me on how to convert shape files into .asc files? I want to convert layers from shape files into this format but I don't know how to go about it in qgis.


texts: Tables within the Reclassify tool- what is it, where is it?

texts:
ArcMap 10.1 here. I have a raster set I'm looking to reclassify. The field I need to re-classify is currently a string of words like "Thing 1" and "Thing 2" and I'm going to set them to be 1 and 2.

I understand that you can create essentially a table within the reclassify tool, save it, and then use it again in another operation. What type of file is this, and where is it found? I can't find it in windows explorer, but it does show up in arc catalog.

Reason is, I have an excel sheet where all of my "Thing 1"s are matched up with the number I want. Let's assume there are thousands of "Thing 1"s and renumbering them within the reclass window would be impractical.

Is there a way to
(A) edit a reclass table within excel
(b) have arcmap reclass tool read in my excel, dbf, csv, or some other kind of file?

Thanks for thoughts, suggestions

D


texts:
My goal:

Store the attributes of a geojson in a csv file
Allow a non-GIS user to edit a cloud-hosted csv file which then updates a webmap
Have the web map be a static app with no GIS/PostGIS/SQL etc server involved

My method:

Load a csv (from say dropbox) and a geojson stored on the webserver
Join the csv to the geojson file using javascript
Style the polygons using the field added from the csv

This works, but it is too slow...
The loading of the csv happens in about 1 second, and the geojson tiles get created almost instantly. But the join function using eachLayer() takes between 7 and 10 seconds - on a desktop machine with a fast internet connection. On a phone performance for the join is awful, maybe 30-50 seconds.
Is there a faster/better way to join the geojson to the csv?
I could be going about this in the absolute wrong way altogether.
I am using leaflet 0.7.7. The geojson and csv contain 4500 features/rows.
Using performance.now() after the Papa.parse, eachLayer and after the tile script gives me this:
csv loaded in 889.39
csv joined in 7406.235000000001
geojson tiles created in 7523.535000000002

And without the join:
csv loaded in 918.0150000000002
geojson tiles created in 1030.2800000000002

Here is the code (wellZoneParcelsData is the omnivore variable and wellZoneParcels is the L.geoJson):
wellZoneParcelsData.on('ready', function () {
  Papa.parse('assets/data/wellston_zone_codes.csv', {
    download: true,
    header: true,
    complete: function(results) {
      var now = performance.now();
      csvloaded = now - start;
      console.log('csv loaded in ' + csvloaded);
      wellZoneParcels.eachLayer(function(layer) {
        var pid = layer.feature.properties.PARCELID;
        var csvData = results.data.filter(function(data){return data.PARCELID == this;}, pid);
        layer.feature.properties.zonecode = csvData[0].zonecode;
      });//end eachlayer function
      now = performance.now();
      csvjoined = now - start;
      console.log('csv joined in ' + csvjoined);
    }//end complete function
  });//end papa parse
});


texts: I am trying to combine L.MarkerClusterGroup with L.TimeDimension: basically, I would like to load a geoJSON and control the appearance and disappearance of points with L.TimeDimension, but in such a way that the points do not appear as individual points but are clustered with L.MarkerClusterGroup or using a different method if simpler. Have you ever seen anything like that? Can you give me some advice or refer me to some example?


texts: Selecting coastal census tracts in ArcMap?

texts: What tool should I use to select or identify geographic units that are adjacent to the U.S. coastline feature in ArcMap?


texts:
I have a coordinate (longitude, latitude) which I get from my map (I'm using OSM with OpenLayers) so I get something like this lon: -46.502541777468984 lat: -23.474943422638063.

What I need? I need to know if there is some function in PostGIS that I can use to search on my database which will return the nearest result that I have in DB.

I know I can do this using a select, but it take too much time.


texts: Connect to remote Sqlserver using arcobject

texts:
I want to connect to a database assigned to a vps with an IP address. I used this string to connect to Database when it was on localhost. How I change it to connect to remote one?

IPropertySet propertySet = new PropertySetClass();
propertySet.SetProperty("INSTANCE", "sde:sqlserver:CRAZYVIRUS-HP");
propertySet.SetProperty("DATABASE", "land_gis");
propertySet.SetProperty("USER", "LandManagment");
propertySet.SetProperty("PASSWORD", "pass");
//propertySet.SetProperty("VERSION", "dbo.DEFAULT");
//propertySet.SetProperty("authentication_mode", "OSA");

workspace = workspaceFactory2.Open(propertySet, 0);

I use ArcObjects 10.2.1 version and C#


texts: Frequency, summarize tool options of ArcGIS for Desktop?

texts:
I'm looking for some general frequency, summarize tool.
What I basically want is some tool like the summarize tool in ArcGIS (right click in the attribute table> summarize) that runs on all attributes and objects that I select (it would be great if that tool could output to excel for example but that is totally optional).
I have way to many objects and attributes to use do this manually...

I know this should be possible to create with ArcPy (although I don't find summarize in Arcpy), but this looks like something that can be useful to a lot of people so I believe something like this must already exist? ArcGIS is no requirement, if it runs on shape files or anything I can export to, that's fine.

If this does not exist, then I will try to develop it myself and maybe post it somewhere online cause it looks like something useful.

EDIT: I've written a small tool in python to basically do summarize. It works on shape files and summarizes to Excel (but other formats are easy to use too) and should be OS independent. If anyone would like to see it I will clean the code up and publish it.

EDIT2: The tool can be found here: https://bitbucket.org/tcoopman/shapesummarize/


texts: Permisson denied when running script or Error 4: no such file or directory(standalone)

texts:
I started working with QGIS, because I need to build a forest fire risk map for a university project. I tried to build a little standalone application to calculate the NDVI. I'm following the python cookbook and made a mix of the standalone app code of chapter 1 and the code to calculate de NDVI. Unfortunately, I get an error I can't find a solution to. I guess it’s because the file of the raster calculator is not being generated properly.
So, when I try to access the layer afterwards I get an error because the file wasn't created.
I have tried a lot of different approaches but I haven't found a solution so far.

I use a Windows PC and Qgis 2.18 Las Palmas.

from qgis.core import *
from PyQt4.QtGui import *
from PyQt4.QtCore import *
from qgis.analysis import *
import sys

app = QgsApplication([], True)
app.setPrefixPath("D:/OSGeo4W64/apps/qgis",True)
app.initQgis()
rasterName="Banda"
raster = QgsRasterLayer("C:/Users/Name/workspace/MyMap/Datos/DatosSatelites/Landsat8_B4.TIF", rasterName)

if not(raster.isValid()):
    QMessageBox.information(None, "DEBUG:", "Error al cargar")

ir = QgsRasterCalculatorEntry()

r = QgsRasterCalculatorEntry()

ir.raster = raster

r.raster = raster

ir.bandNumber = 2

r.bandNumber = 1

ir.ref = rasterName + "@2"

r.ref = rasterName + "@1"

references = (ir.ref, r.ref, ir.ref, r.ref)

exp = "1.0 * (%s - %s) / 1.0 + (%s + %s)" % references

output = 'C:/Users/Name/workspace/MyMap/Datos/Output/ndvi1.tif'

e = raster.extent()

w = raster.width()

h = raster.height()

entries = [ir,r]

**ndvi = QgsRasterCalculator(exp, output, "GTiff", e, w, h, entries)**

ndvi.processCalculation()

**lyr = QgsRasterLayer(output, "NDVI")** #in this line I get the message

algorithm = QgsContrastEnhancement.StretchToMinimumMaximum

limits = QgsRaster.ContrastEnhancementMinMax

lyr.setContrastEnhancement(algorithm, limits)

s = QgsRasterShader()

c = QgsColorRampShader()

c.setColorRampType(QgsColorRampShader.INTERPOLATED)

i = []

qri = QgsColorRampShader.ColorRampItem

i.append(qri(0, QColor(0,0,0,0), 'NODATA'))

i.append(qri(214, QColor(120,69,25,255), 'Lowest Biomass'))

i.append(qri(236, QColor(255,178,74,255), 'Lower Biomass'))

i.append(qri(258, QColor(255,237,166,255), 'Low Biomass'))

i.append(qri(280, QColor(173,232,94,255), 'Moderate Biomass'))

i.append(qri(303, QColor(135,181,64,255), 'High Biomass'))

i.append(qri(325, QColor(3,156,0,255), 'Higher Biomass'))

i.append(qri(400, QColor(1,100,0,255), 'Highest Biomass'))

c.setColorRampItemList(i)

s.setRasterShaderFunction(c)

ps = QgsSingleBandPseudoColorRenderer(lyr.dataProvider(), 1,s)

lyr.setRenderer(ps)

QgsMapLayerRegistry.instance().addMapLayer(lyr)

exitcode = app.exec_()

QgsApplication.exitQgis()

sys.exit(exitcode)


texts:
I have problem with getting file path outside the class, my method does not work and I don't know why.
class ButtonImport(QPushButton):
    def __init__(self, accepted_formats):
        super().__init__("Import")
        self.__accepted_formats = accepted_formats
        self.__filepath = None
        self.clicked.connect(self.__handle_select_file)

    def get_filepath(self):
        return self.__filepath

    def __handle_select_file(self):
        self.__filepath, _ = QFileDialog.getOpenFileName(self, "Select file")
        file_extension = path.splitext(self.__filepath)[1].lower()
        if file_extension not in self.__accepted_formats:
            print("Please choose file one more time")


texts: Leaflet Map Clustering zoom in not working on close points

texts:
The Cluster don't get zoom in, it just split like this if the points are close to each other , in other clusters if the points are not closer they just get zoom and working fine, I want to disable this affect if possible and make it zoom on clicking on it.Kindly guide me if this is possible or not.

Points :

LAT 31.4926473    LON 74.3759229
LAT 31.4923959    LON 74.375927

   function setMap(cat_ids) {

    $.ajax({
        url: '{!!route('all-outlets-map-client-portal')!!}',
        data: {
            cat_ids: cat_ids, project_filter_applied:project_filter_applied, geofilters:geofilters
        },
        dataType: 'json',
        success: function(data) {
            locations = data;
            ids       = data.project_ids;
            // console.log(ids);
             initialise = 1;
            populate(cat_ids);
        }
    });

    function populate(cat_ids) {

        $.each(cat_ids, function( index, value ){
            cities = '';
            layers_variables[value]           = L.layerGroup();
            layers_variables_clusters[value]  = new L.MarkerClusterGroup();

            map.addLayer(layers_variables_clusters[value]);
            var color = getRandomColor();
            for (var i = 0; i < locations[value].length; i++) {
                var id  = locations[value][i].id;
                var lat = locations[value][i].latitude;
                var lng = locations[value][i].longitude;

                var marker = new L.circleMarker(new L.LatLng(lat, lng), {
                    radius: 6,
                    fillColor: color,
                    color: color,
                    fillOpacity: 1,
                    renderer: myRenderer
                }).addTo(layers_variables[value]);

                marker.on('click', onClick_Marker)
                // Attach the corresponding JSON data to your marker:
                marker.id     = id;
                marker.cat_id = value;
                marker.color  = color;
                markerList.push(marker);
                layers_variables_clusters[value].addLayer(marker);
            } //end for loop
        });
    } //end Populate function
    if(initialise == 0){
        L.control.layers(baseLayers).addTo(map);
    }
}


texts:
I have an old Mapnik's XML file, from which I need to serve out data as a WMS service.

I've gone through the Mapnik site, and I didn't find anything. I've also looked at Mapnik's OGC server, and looks like there has been no development in the last 10 months.

How do I serve out a WMS using mapnik?


texts: ArcObjects SDK for .NET JSONToFeatures fails when large JSON file is used

texts:
I have a large JSON file which I download from ArcGIS Server REST (size is 140MB). When I try to execute the JSONToFeatures Geoprocessing tool on this file I get an error.

Code: (call below method with arguments toolName="JSONToFeatures", a path for input json and a path for output SHP file).

    private object RunTool(string toolName, object[] parameters) {
        IVariantArray gpParams = new VarArrayClass();
        IGeoProcessorResult gpResult = new GeoProcessorResult();
        foreach (object o in parameters)
            gpParams.Add(o);
        object rev = new object();
        try {
            gpResult = Geoprocessor.Execute(toolName, gpParams, null);
            return gpResult.ReturnValue;
        }
        catch (COMException comEx) {
            throw new Exception(Geoprocessor.GetMessages(ref rev), comEx);
        }
        finally {
            ReleaseObjects(new object[] { gpParams, gpResult });
        }
    }

When executed returns the following error:
Executing: JSONToFeatures input_path output_path
Start Time: Thu Jan 31 14:32:43 2019
ERROR 999998: Unexpected Error.
Failed to execute (JSONToFeatures).
Failed at Thu Jan 31 14:32:45 2019 (Elapsed Time: 1,48 seconds)

When I execute JSONToFeatures tool from within ArcDesktop ModelBuilder execution proceeds successfully.

Please note that on small files of size 1-10 MB (which is the normal size of files processed by the code above) there are no problems.

This error is very generic. How do I need to proceed?


texts: Is there a programatic method to simplify a polygon?

texts:
I have an array, containing lat/lon coordinates which form a polygon(the boundaries of certain cities on the map). I got this data from OpenStreetMap's database, and my problem is what follows:

Certain location boundaries, are formed by 5-6000+ lat/lon coordinates, which is a very large data for me to work with. In the meanwhile, I've observed that a lot of these polygons could be trimmed down to contain a fewer points, because a lot of lat/lon coordinates are placed on the same line.

What I would be interested in, is if it would be possible, to get a polygon, and programatically check if multiple points of it reside on the same line(with a small error margin), and delete unneeded(in my case) points, thus so minifying the polygon itself?


texts: Coordinates of upper left corner in the world file?

texts:
According to Wikipedia and ESRI last two lines in the world file represent x and y coordinates of center of upper left pixel.

Now, if I have a world file like this:

32.0
0.0
0.0
-32.0
691200.0
4576000.0

...can I say that the "real" upper left corner of the image (upper left corner of the upper left pixel) is:

X = 691200.0 - (32.0 / 2) = 691184.0
Y = 4576000.0 - (-32 / 2) = 4576016.0

...and should I use these coordinates if I want to use the image with Bing Maps or Google Maps (after converting to Lat/Long off course), or should I stick with original ones (pixel center coordinates).

UPDATE:

I just realized that the .tab file associated with the map has coordinates of the upper left pixel (0,0) calculated the same way as in my question above. Not sure about the exact formula but resulting coordinates are the same - shifted left for half of the pixel width and up for half of the pixel height. Now I'm pretty confident that this is the right way to do it.


texts: Programatically setting default Geodatabase and Workspaces for map document?

texts:
I have a QC process that has a bunch of models that need to be sequentially run, and then the results manually inspected between each model.

I know how to manually set the default GDB for the Map Document (File-->Map Document Properties-->Default Geodatabase).

I also know how to manually set the default Geoprocessing Workspaces (Geoprocessing-->Environments-->Workspace).

I already have a "setup script" that generate a GDB, Feature Datasets, define the spatial reference, etc, and I would like to be able to integrate these steps into the setup procedure as well.

Can these settings be automatically set?  For example, I have tried arcpy.env.workspace, but that only appears to work within a script, not setting it at the "application" level.  Is this possible?

I know it seems like a small thing, but I have many "areas of interest" that we have to look at, and each one gets placed in its own MXD and GDB (so yes, I could do it by hand... over & over again).


texts: How to intersect shapefiles with elevation models?

texts:
I have a high resolution digital elevation model (DEM) and a line shapefile displaying a bicycle track. The DEM has 1 meter resolution. The line shape is something like 25 km long.

What I need is a detailed elevation profile of that track, let's say one coordinate every 1 meter.

What I've tried so far is to convert the shp file to gpx files and used tools like gpsvisualizer or gpxplot to plot elevation profiles of the track. I used a great circle calculation afterwards to find out the distances between the points.

The distance between single points varies from 3 to 700 meters which makes this plot useless. The main problem is that the conversion of line shape to GPS coordinates only generates the corners of the line but not points in-between somewhere on the edge of the line.

I'm wondering is there an easy way how to use my elevation model (available in AutoCAD files format adf) to intersect it somehow with the line shape to plot some detailed elevation profile both in numerical (lat, lon & elev. or x, y & z) and graphical ways?

I need this to do further analysis on the given bicycle track.


texts: Consolidating the attribute data from underlying smaller polygons to a bigger layer on top

texts:
I am using ESRI's ArcMap 10.2

I have an underlying polygon layer made up of small pieces and they have census data attached to them (think of counties for the state of Texas).

I want this data averaged and visible in the attribute table for the larger overlying polygon (think of the state of Texas), therefore I want to have the average census data across these counties visible as one row of data for Texas. E.g, Average income from the overlying counties, averaged out for the entire overlying polygon.

Kinda looks like this:

Polygon - AverageIncome - WhiteRace - PeoplePerHousehold

Polygon1 - 24124        - 13410     - 2.6

Polygon1 - 29481        - 25281     - 3.2

Polygon1 - 68128        - 43591     - 4.2

And I want it to end up like this for the polygon that covers them:

Polygon - AverageIncome - WhiteRace - PeoplePerHousehold

Polygon1 - 40577.6        - 27427     - 3.3

For the life of me I have spent hours searching the internet and trying different tools, without a solution.


texts: Skip duplicate records when appending using arcpy?

texts: The objects I want to check have unique ID's (edit: these ID's are different from OID's and will remain the same when copied). I'm appending one file geodatabase to another, and I want to make sure I'm skipping any duplicate objects. What would be the best way to go about this?


texts: Correctly Transforming (Unprojecting?) from Lambert Conformal Conic?

texts:
I started out by posting https://stackoverflow.com/questions/6963991/how-can-i-transform-the-coordinates-of-a-shapefile

In short, I found a shapefile here, and I want to be able to convert that to Latitude/Longitude (Geodetic?) coordinates that I can store in my PostGIS database.  I am having trouble working with the Ruby libraries, but that is a separate issue.  First, I need to correctly understand the projection.

In the linked zip file, there is a .prj file, which contains the following (It is only one line, but I broke it up for readability):

PROJCS["NAD_1983_StatePlane_California_III_FIPS_0403_Feet",
GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",
SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],
PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",6561666.666666666],PARAMETER["False_Northing",1640416.666666667],
PARAMETER["Central_Meridian",-120.5],PARAMETER["Standard_Parallel_1",37.06666666666667],
PARAMETER["Standard_Parallel_2",38.43333333333333],PARAMETER["Latitude_Of_Origin",36.5],
UNIT["Foot_US",0.3048006096012192]]

I was unable to find tools that translate .prj files into a projection string, so I used this as a reference and came up with:

+proj=lcc +datum=NAD83 +ellps=GRS80 +lat_1=37.06666666666667 +lat_2=38.43333333333333 +lat_0=36.5 +lon_0=-120.5 +x_0=6561666.666666666 +y_0=1640416.666666667 +units=us-ft +pm=greenwich +no_defs

I happen to know that the points in this shapefile should correspond to locations in San Francisco, so I could test my work:

cs2cs -f "%.8f" +proj=lcc +datum=NAD83 +ellps=GRS80 +lat_1=37.06666666666667 +lat_2=38.43333333333333 +lat_0=36.5 +lon_0=-120.5 +x_0=6561666.666666666 +y_0=1640416.666666667 +units=us-ft +pm=greenwich +no_defs +to +proj=lonlat +datum=WGS84 +ellps=WGS84
6011287.4999795845 2100857.2499904726
-164.18117482   17.72087079 0.00000000

That's about 1000km west of Hawaii, so clearly that's not right.  I tried many different things, but the closest I got was by stripping the +units=us-ft off.

cs2cs -f "%.8f" +proj=lcc +datum=NAD83 +ellps=GRS80 +lat_1=37.06666666666667 +lat_2=38.43333333333333 +lat_0=36.5 +lon_0=-120.5 +x_0=6561666.666666666 +y_0=1640416.666666667 +pm=greenwich +no_defs +to +proj=lonlat +datum=WGS84 +ellps=WGS84
6011287.4999795845 2100857.2499904726
-126.98864867   40.47509388 0.00000000

That's about 500 km NW of San Francisco.

What am I missing here?
If I import the shapefile into the online ArcGIS Explorer, it understands it perfectly.


texts: How to get the best accuracy of coordinates for web GIS- application?

texts:
I would like to get advice from you about one thing. I would like to develop my own web GIS application. It would be about some sport places in my hometown; it would be good practice for my education - I'm still a student. First, I need to get coordinates of the features, then I think to put it in AutoCad, then export to .shp file, then work with it in QGIS and use Postgres and PostGis for databases and than use Mapserver and so on.

But I need help with the first step, getting the coordinates. I don't have any high accurate instrument like some good GPS, so is there any chance to get coordinates of the places by Google Earth? How accurate is that and is it even possible to do it like that? I'm looking forward to someone's answer and please don't get mad if I said something stupid because as I said, I'm in the beginning of getting to know all this GIS web development and things like that. Thanks!

Thank you all for your answers, it all helped me in some way. The only thing that I'm not sure is accuracy of osm and bing maps. I was editing some stuff in JOSM recently and I've realized that bing maps don't really fits in the right place if you compare it to gpx data. (there is at least 2m difference) I'm not still quite sure is it accurate enough for me. I would like to achieve 1m maximum level of accuracy.. :/


texts: Visual Studio 2015 Error when running ArcMap add-in

texts:
I'm running a simple ArcMap add-in in Visual Studio and am getting this error:

  Severity  Code    Description Project File    Line    Suppression State
  Error     The "ValidateAddInXMLTask" task failed unexpectedly.
  System.IO.FileNotFoundException: Could not load file or assembly
  'Microsoft.VisualStudio.Shell.12.0, Version=12.0.0.0, Culture=neutral,
  PublicKeyToken=b03f5f7f11d50a3a' or one of its dependencies. The
  system cannot find the file specified. File name:
  'Microsoft.VisualStudio.Shell.12.0, Version=12.0.0.0, Culture=neutral,
  PublicKeyToken=b03f5f7f11d50a3a'    at
  ESRI.ArcGIS.AddIns.SDK.MSBuild.ValidateAddInXMLTask.Execute()    at
  Microsoft.Build.BackEnd.TaskExecutionHost.Microsoft.Build.BackEnd.ITaskExecutionHost.Execute()
  at
  Microsoft.Build.BackEnd.TaskBuilder.d__26.MoveNext()

  WRN: Assembly binding logging is turned OFF. To enable assembly bind
  failure logging, set the registry value
  [HKLM\Software\Microsoft\Fusion!EnableLog] (DWORD) to 1. Note: There
  is some performance penalty associated with assembly bind failure
  logging. To turn this feature off, remove the registry value
  [HKLM\Software\Microsoft\Fusion!EnableLog].   FirstAddIn

I'm using ArcGIS for Desktop and ArcObjects SDK for .NET version 10.4. I obtained my Arc license and SDK via ISO on MyEsri and for this add-in I'm using .NET 4.5.

My issue seems very similar to this one as I also got the same dialog box pictured there, but in the accepted answer it has me check Use Managed Compatibility Mode and it was grayed out when I tried it. I also had no luck with the other answers. I'm new to Visual Studio.


texts:
I am working with the vector files:

Landuse classification shapefile for a certain area with 4 landuse categories - builtup, vegetation, water and open spaces. (area bounds = 95sqkm). I classified a satellite image of Sentinel-2 in Google Earth Engine and converted the raster classified file to vector in ArcMap.
Grids against the same area bounds (grid size 50m x 50m)

I want to get the grids that have an overlapping of 50% or more area with the Built-up landuse from the landuse file.

Is there any tool in ArcMap that can do the job efficiently and quickly?


texts: How do you validate longitude and latitude coordinates in Python?

texts:
I am working on a bunch of shapefile data, and I need to process it and after some kind of processing, I need to feed the data to an indexer (in this case ElasticSearch). The problem is I just discovered that there are some anomalies in the data. I found this in one of my polygons.

[(593.0587471450001, -1.8046608309999783),
 (593.0587471160001, -1.8046608379999611),
 (593.0587470810001, -1.804660872999932),
....]

Now correct me if I am wrong, but aren't longitudes suposed to be in the range of -180 to 180 and latitudes in the range of -90 to 90? So the numbers above just don't make sense. Is there any function in shapely, GeoPandas or any other library that can validate/check whether a coordinate is valid or not? Or should I do that manually?

UPDATE:

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

Also I don't know if this is relevant, but I only found this one polygon with very big numbers in it. The other polygons in the same file are fine (in the right range).


texts: Layers are not visible after using a custom orthographic projection in QGIS 3.10

texts:
I was planning to make a global orthographic map, that would show country borders and oceans. Therefore,i followed this method which plots a map similar to my requirements. However, i modified the codes so that the map is centered on Asia:
+proj=ortho +lat_0=10 +lon_0=90 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs

Although the land areas were projected just fine, i was having problem with the ocean's layer, as you can see in these pictures:
I decided to change the projection of the ocean layer to the custom projection, but that makes the ocean layer entirely invisible:

There has been a similar post regarding this problem, where it was suggested to use uncheck the simplify geometry options in the layer rendering, but that didn't work for me. I'd be really grateful if someone suggested me any way to properly project the ocean layer in the globe, like in this picture:  Thanks in advance.


texts:
When I use the stratigraphy plot function in Midvatten on an East Anglian coast quaternary data set (West 1980) that I have been digitizing, I find that the upper layers of strata disappear at the top of the window as I expand the window. Also not all the strata get the right colour coding and some strata appear out of sequence. Here is a screen shot example with relevant attribute tables:

Hopefully the posted screen shot will be expandable so the attribute tables are legible. Chalk should be the base layer in borehole HC but it comes out as the second layer in the plot and only 7 out of 13 strata are visible.

I have:

a) checked that the window chopping off the top of the plot is not associated with the hi-resolution screen by trying the same plot on a computer with a lower resolution screen.
b) renumbered 'stratid' to start every 'obs_point' from 1 (originally I numbered the 'stratid' consecutively for the whole dataset 1-1370).
c) cut the number of facies from 74 geoshort to 25 strata with the hope that simplification might help.(Josef Kallgarden  suggested a and b).

What else might I try to fix these issues?


texts:
I am looking at a way to specify a different nodata value for each band in rasterio.

I tried to:

use the nodata key in the profile with a tuple or list of nodata values,
set the nodatavals attribute,

but none of these methods worked!

Is that feasible?


texts:
Ive updated my spatial_ref_sys table in postgres with spatial reference statement from http://spatialreference.org/ref/esri/102737/.

I can perform the tranformation from srid 4326 to srid 102737 and receive no errors...

INSERT into spatial_ref_sys (srid, auth_name, ...

SELECT id, name, ST_Transform(subs.geom, 102737) as geom
INTO subs4201

However the results when viewed on a map are displaced.

I can test this by using arcgis to reproject the data and the results overlay the wgs84 original locations....

How can I correct for this?

Edit: The displacement is approx. 1/2 mi.

Edit2: Adding ArcMap image as requested by @mkennedy


texts: PostGIS + pgRouting: Adding dummy links to a road network and performing routing analysis

texts:
I'm a beginner when it comes to postgis and pgrouting development.

For a project I'm working on and a bit of the algorithm:

I have a  polygon layer that represents zones
I have a line layer that represents roads based on classifications: Freeway, Arterial, and Collector, and Local
I created a point layer of all the zone centroids
I would like to create a table of  n x "dummy link(s)". These links connect the centroid to the road layer based on the road classfication (above), the link is within the polygon, and based on nearest neighbour (nearest proper road classification segment to centroid).
A dummy link has infinite capacity and has a cost of zero.
These links should be added to their own table
I would like to start the routing from the centroid layer.
At each of the centroids, I would like to do a minimum spanning tree or Dijkstra to all other centroids of the network.
The end goal would be a matrix of some sort.
I'd like this operation to be somewhat fast, as I would have about 200 zones (400 once you consider two way traffic), with a road network of roughly 1,000,000 links. Road weights/costs will change as they are time of day dependent; therefore, there will be at least 24 iterations per zone.
If possible, I'd like to limit most of my development to python as I'm comfortable with python and libgdal.

Question:

Is this feasible using postgis + pgrouting, especially considering that I am not proficient at SQL?
Can you critique my model, suggest improvements, "gotchas", anything that would make my life a little bit easier?


texts: ENVI files (datafile +hdrfile) not displaying properly in python

texts:
I am trying to topographically correct my already atmospherically and radiometrically corrected images. To do so I need to calculate illumination conditions at the time of image acquisition. Terms in this equation are slope and aspect. The two files are stored in a raster ENVI file format, generated from a ASTER-II DEM in ArcGIS. There are other lines beforehand where i define wdir and other variables, but assume that everything works super smooth until this point. Code is:

rgb = np.stack([refl_surf[5,:,:],refl_surf[3,:,:],refl_surf[2,:,:]],  axis=2)
imshow(rgb)
show()

def enviread(datafile, hdrfile):
   info = envi.open(wdir+hdrfile, wdir+datafile)
   data = info.load()
   return data, info

[slope,info1]=enviread('Slope_DEM_UT11.dat','Slope_DEM_UT11.hdr')
[aspect,info1]=enviread('Aspect_DEM_U11.dat','Aspect_DEM_U11.hdr')

fig = plt.figure()

plt.subplot(121)
imshow(np.squeeze(slope))
colorbar()

plt.subplot(122)
imshow(np.squeeze(aspect))
colorbar()

plt.show()

When showing the images, instead of some topography I see just the following yellow image:

This is wrong because when I do it with the slope and aspect maps, provided by my teacher, running the same code they display like so:

It must be something I've done while generating the slope and aspect but what can it be?

Slope and aspect display without any problems in ArcGIS.

DEM Coordinate system was converted from its original geographic one to the UTM zone where my study area is before extracting any of the derivatives.

It turns out I was looking not at the exported ENVI files, but at my original images. Apparently exporting to an ENVI file format messes the images up, they show just a black monotone image even in ArcGIS.

Is this some kind of inherent flaw of exporting to ENVI?


texts: Visually Merging ZIP polygons and Extracting by boundary from SQL Server 2008

texts:
i have a lot of shapefiles. For simplicity, lets say all the zipcodes of the USA. I've imported them into SQL SERVER 2008.

I wish to merge lots of the shapes into one shape.

I'm sorta struggling to do this with Sql Server 2008 and sql script.

Are there any programs out there that can do this? eg. click on these shapes and merge into one?

Also, can i cut shapes out? So lets imagine i have a whole pile of zipcodes together .. and then i draw my own shape .. can I cut out my shape from those zipcode leaving a nice whole.

Once more I can sorta do this with SQL Script but I'm getting really tired of this and would really love it if this stuff is all visual.

Any ideas, folks?


texts: Geotools 8.4 CRS.decode invalid(?) reference to HTTP_URI_OGC

texts:
I'm trying to do a simple CRS.decode to lookup a CoordinateReference system:

E.g. CoordinateReferenceSystem target_crs = CRS.decode("EPSG:4326");

However this fails with a reference to the Citations.HTTP_URI_OGC field, and an invalid field error. According to the documentation this field is introduced in v9?

Anyone else experienced this, or got an idea of what I'm doing wrong?

Regards

Rune Sørensen


texts:
arcpy: I would like to join a .xls table to a .shp table using two common fields, instead of one common field.

Is this possible in ArcGIS using python?


texts: Hec-Ras 5.0: Adding aerial imagery to geometric view

texts: I'm having trouble displaying an aerial imager in the HEC-RAS geometric view in version 5.0.1. I've added the aerial, a project TIF file, to the RAS Mapper and it shows up without issue. When I go to the geometric view and select "select layers to view in background" and select the imagery, nothing shows up. I'm in a blank geometry file without any data created yet. I've tried doing View > Set Geometry Extents, but that still leaves with with white space.


texts:
I am looping through Municipalities and clipping a layer called AllZones containing Multiple Facility (FACID) with multiple chemicals per Facility (CID). 1 to many relationship

           FacilityID,        CID,  other rows......
ex rows:        001              2
                001              2
                001              2
                001              5
                001              5
                001              7
                001              7
                002              11
                002              12
                002              14

that result layer I then dissolve it down just by the FACID (call it labeldissolve) so I can generate 1 label for each facility on the map based on the new objected id.

   OBJECTID,  FacilityID
       1          001
       2          002

My next step is to create a report of the AllZones layer info for that given municipality. I want those created OBJECTID labels in the labeldissolve layer to correspond to the report. I tried joining the two tables in ArcMap based on the FACID and it works fine.

When I tried it in arcpy I did not have the option to choose keep only matching records

here is my code

arcpy.JoinField_management("AllZones", "FACID", "LabelDissolve", "FACID","")

should I be using a spatial join?


texts: Shapefile is saved in EPSG:4326 despite being set to EPSG:32651 via geopandas

texts:
This is the code

def loadPoint():
    df = pd.read_csv(r'C:\Users\imper\Documents\new_GISfiles\earthquakesPH\emscPhilippines2008to2020_lonLat.csv')
    geometry = [Point(xy) for xy in zip(df.Longitude, df.Latitude)]
    crs = 'GEOGCS["WGS 84",' \
          'DATUM["WGS_1984",' \
          'SPHEROID["WGS 84",6378137,298.257223563,' \
          'AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]], ' \
          'PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],' \
          'UNIT["degree",0.0174532925199433,' \
          'AUTHORITY["EPSG","9122"]],' \
          'AUTHORITY["EPSG","4326"]]' \
          'PROJECTION["Transverse_Mercator"], ' \
          'PARAMETER["latitude_of_origin",0],' \
          'PARAMETER["central_meridian",123],' \
          'PARAMETER["scale_factor",0.9996],' \
          'PARAMETER["false_easting",500000],' \
          'PARAMETER["false_northing",0],' \
          'UNIT["metre",1,' \
          'AUTHORITY["EPSG","9001"]],' \
          'AXIS["Easting",EAST],' \
          'AXIS["Northing",NORTH],' \
          'AUTHORITY["EPSG","32651"]]'
    geo_df = GeoDataFrame(df, crs=crs, geometry=geometry)
    geo_df.to_file(driver='ESRI Shapefile', filename='eqpts.shp')

def plotPoint():
    basemap = gpd.read_file(r'C:\Users\imper\Documents\new_GISfiles\earthquakesPH\PH_provs.shp')
    pt = gpd.read_file(r'C:\Users\imper\Documents\new_GISfiles\earthquakesPH\eqpts.shp')
    ax = basemap.plot(cmap = 'jet', column = 'NAME_1', figsize=(10,10))
    pt.plot(ax=ax, marker="o", color="red", markersize=5)

loadPoint()
plotPoint()

But CRS is wrong

pt.crs
Out[44]: {'init': 'epsg:4326'}

CSV is from EMSC (earthquaks location)

  Latitude,Longitude

  13.81,121.1

  13.76,120.92

  10.45,125.2

  6.69,125.23


texts: Error "Element.get: Parameter 'object' is required" mapping to a list in Google Earth Engine

texts:
I'm running my code where I create weekly averages of a variables which are mapped to a list of days of year (DOY).
var SSTA_7d =  ee.ImageCollection.fromImages(
  years.map(function (y) {
      return DOY.map(function(d) {
      d = ee.Number(d);
      var SSTA_date = series_daily.filter(ee.Filter.calendarRange(y, y, 'year'))
                      .filter(ee.Filter.calendarRange(d, d.add(7), 'day_of_year'));

      var img_last = ee.Image(SSTA_date.sort('system:time_start',true).first());
      var size = SSTA_date.size();

      var SSTA_filtered_7d = ee.Image(ee.Algorithms.If(size.eq(0),
                                            ee.Image.constant(0).toFloat(),
                                            ee.Image(series_daily.filter(ee.Filter.calendarRange(y, y, 'year'))
                                            .filter(ee.Filter.calendarRange(d, d.add(7), 'day_of_year'))
                                            .mean()
                                            )));

      return SSTA_filtered_7d.set('year', y)
              .set('DOY', d)
              .set('system:time_start', img_last.get('system:time_start'));

  });
  }).flatten()
);

But I'm getting this error.
ImageCollection (Error)
Element.get: Parameter 'object' is required.

From reading it might be because the image collection I'm using is empty or a server side issue. But I can't quite figure out where the issue lies. or why the image collection is empty if that's the case. My full code is here. Lines 118-141.


texts:
In ArcToolbox I have a script (arcpy) that calculates some info plus a date into a database table.  The date required to write is often (but not always) the current date, so the tool asks for input of the date to write.

I'd like to have the Tool automatically populate the current date into this input field so that the user only needs to modify the date on the occasions where the date is not today.

How can I set the Default Value to be the current date?  I only seem to be able to enter dates into the Default field in the tool properties (as it's a Date field), so don't seem to be able to use anything like date.today() or similar.


texts: We writing pixel values in Google Earth Engine

texts:
I am working with the Ocean Color SMI: Standard Mapped Image MODIS Aqua Data dataset.
Using the Sea Surface Temperature band, I am trying to assign the value "0" to any pixel in my region that has a value less than 30. Conversely, I need to assign a "1" to a pixel that has a value greater than or equal to 30. Essentially, some sort of binary process on the pixel data so that then when I visualize the map, there will be only two colors, the "0" and the "1."
How can I achieve this?


texts: Hosted feature layer not displaying in Scene Viewer - ArcGIS Online

texts:
I have a scene in ArcGIS Online that I had managed to display a hosted feature layer of polygons on but today I messed up my contents and now the feature layer will no longer display in the scene. I didn't change anything about the feature layer itself (I deleted a few unused items to free up storage, one of which was a tile layer made from the feature layer, but I checked in advance and nothing I deleted was linked to the scene that I could tell) I recreated a tile layer and a vector tile layer but neither of them work in scene, I just get an error message so that can't be it.
I can add other polygon hosted feature layers that appear to be configured and shared the exact same as the one I want, they show up automatically in the default style 2D 'on the ground'.
Anything I can do to make the one I want show up? I remember when I was learning how to add it I had trouble and then an 'aha!' moment but I can't find any help online and can't figure out how to make it work this time.


texts:
We are working on our new vehicle tracking system. We have been making use of mySQL for the past two years. However we have been discussing moving to a noSQL database.

We have a high number of reads and writes. For reads we have been told that mongodb will take up to much of our available ram on our servers. So which noSQL will best suit our needs, without causing us any furthur problems? Or is it better to hang on to mySQL?

Thanks in advance!


texts:
I am not a GIS professional, only a casual orienteering enthusiast seeking an explanation from experts.

The plat of my property gives bearings and distances for each of the lines, as does the deed, however the bearings are wildly off from either true north or magnetic north.  It specifies that "course bearings are based on Connecticut Coordinate System (NAD 27)", and also displays an arrow showing the direction of grid north ("Connecticut Coordinate Grid North") about 20 degrees clockwise from a parallel of the edge of the page, a far cry from the degree or so I usually see for grid north.  Also unlike a USGS Topo map, there is no reference to either true or magnetic north, just the line marking grid north.

If I assume that the sheet is oriented to have true north be 'up', then the bearings start to make sense; e.g. an azimuth given as 341 on the sheet ends up being 1-2 degrees as I stand out in the yard and measure it with a sighting compass (correcting for declination).  In checking out the state plane coordinate system for CT, which I admittedly do not profess to fully understand, the correction (theta) for my longitude should only be about 20 minutes, and certainly not 20 degrees.

It's easy enough to just add 20 to each of the bearings, but I am interested in the how and why the surveyors would orient grid north in this way, as well as the mechanism someone would use to determine the exact grid-to-geodetic correction instead of just guesstimating from the edge of the paper.

Is there some missing information that I don't have?


texts: Rotate coordinates from climatic data in raster format (cordex-EU domain 11) to usual coordinate systems

texts:
I am trying to create a dataset that shows future climatic condition in Europe aggregates to NUT2 administration level in R. To do so, i work with cordex data (domain EU 11), that is coming withrotated coordinates.

I used the data hurs_EUR-11_ICHEC-EC-EARTH_rcp85_r1i1p1_SMHI-RCA4_v1_day_20310101-20351231.nc (data here ) which is the relative humidity
Its looks as follows :

 6 variables (excluding dimension variables):
    double height[]
        axis: Z
        long_name: height
        positive: up
        standard_name: height
        units: m
    double time_bnds[bnds,time]
    char rotated_pole[]
        grid_mapping_name: rotated_latitude_longitude
        grid_north_pole_latitude: 39.25
        grid_north_pole_longitude: -162
    float hurs[rlon,rlat,time]
        grid_mapping: rotated_pole
        _FillValue: 1.00000002004088e+20
        missing_value: 1.00000002004088e+20
        standard_name: relative_humidity
        long_name: Near-Surface Relative Humidity
        units: %
        coordinates: lon lat height
        cell_methods: time: mean
    double lon[rlon,rlat]
        standard_name: longitude
        long_name: longitude
        units: degrees_east
    double lat[rlon,rlat]
        standard_name: latitude
        long_name: latitude
        units: degrees_north

 4 dimensions:
    time  Size:1826   *** is unlimited ***
        standard_name: time
        units: days since 1949-12-01 00:00:00
        calendar: standard
        long_name: time
        bounds: time_bnds
        axis: T
    bnds  Size:2

  2 "vobjtovarid4: **** WARNING **** I was asked to get a varid for
  dimension named bnds BUT this dimension HAS NO DIMVAR! Code will
  probably fail at this point"
          rlon  Size:424
              standard_name: grid_longitude
              long_name: longitude in rotated pole grid
              units: degrees
              axis: X
          rlat  Size:412
              standard_name: grid_latitude
              long_name: latitude in rotated pole grid
              units: degrees
              axis: Y

First, i created an array

hurs.array <- ncvar_get(data2, "hurs")

Then i create a grid with the projection as suggested in this post
Choosing CRS / PROJ4 string for EURO-Cordex rotated pole projection?
Note that the issue raised about this post in about an simple feature and not about raster.

r_brick <- brick(hurs.array, xmn=min(lat), xmx=max(lat), ymn=min(lon), ymx=max(lon), crs=CRS("+proj=ob_tran +o_proj=longlat +o_lon_p=-162.0 +o_lat_p=39.25 +lon_0=180 +lat_0=0"))

Now my question is  how to transform this data into the EU projection?
with

r_brickP<- projectRaster(r_brick, crs="+init=epsg:3035")

I get

  not finiteno non-missing arguments to min;

What is going on?

Is this the right approach?

Is there a better approach to unrotate the data?

If yes which one?

What is my most efficient way to then perform a zonal statics with as layer with epsg:3035 projection?

I also tried the function rotate, but i got a line suggesting that this is not this method was not appropriate.


texts: Reserved Ampersand Character in Label Expression of ArcGIS for Desktop?

texts:
I have this code which works...

def FindLabel ([NUMBERNAME], [NUMBER], [BLOCK], [SURNUM]):
      lineone = "A:"+ [NUMBER]
      absName = "<CLR red='255'><FNT size = '14'>" + [NUMBERNAME] + "</FNT></CLR>"
      if [NUMBERNAME] is not None:
          return lineone + '\n' + absName + '\n' + "BLK:" + [BLOCK] + '\n' + "SEC:"+ [SURNUM]

However I have some field values in [NUMBERNAME] that contain the reserved ampersand character (&). If the value contains the ampersand character the label looks like this...

<CLR red='255'><FNT size = '5.5'>T&P RR CO</FNT></CLR>
BLK:
SEC:

Versus (in this example) a label that looks like...

'T&P RR CO'

I have found this documentation, which references this code to fix the ampersand issue....

 Function FindLabel ([LABELFIELD])
  NewString = Replace([LABELFIELD],"&","&amp;")
  FindLabel = "<ITA>" & NewString & "</ITA>"
End Function

This is in VBScript. How would I do something similar to this using the Python parser of ArcGIS 10.3 for Desktop?


texts: How to dock the default ObjectInspector in a new parent container?

texts:
Into my object inspector fever, I noticed that the default Object Inspector object does not dock anymore.

I only have a TabPage control and two tabs, one of those hosting the default inspector. (check black marks)

The other tabs works fine.

Also, this is the code I'm using for the Inspect() method:

    public void Inspect(IEnumRow objects, IEditor Editor)
    {
        SetParent(_Inspector.HWND, DefaultPanelHandle);
        ShowWindow(_Inspector.HWND, SW_SHOW);
        _Inspector.Inspect(objects, Editor);
    }

Why the default does not dock anymore? Any workarounds?


texts:
I would like to be able to convert a single feature (or potentially many features) in an ArcGIS layer into Well Known Text (WKT) format during a geoprocessing operation using Python and ArcPy (no non-ArcPy modules).  The purpose is to then pass the WKT on to SQL Server Spatial and do additional processing outside of the ArcGIS GP toolset.  Is there a method for converting a ArcGIS geometry features to WKT via ArcPy?

I have already read the following, without finding what I am after:

Shapefiles to Text (references a tool in the deprecated "Samples" toolbox)
Convert between ESRI Geometry and WKT (uses ArcObjects or non-ArcPy modules)
How to convert a shapefile to WKT? (Uses QGIS and OGR)

There apparently used to be a tool called "Write Features To Text File" (which appears to be a Python script) that was in the Samples toolbox, but that toolbox was deprecated at version 10 and I cannot find a copy of the Samples Toolbox (I have v10.0) on my machine.  If there are no current solutions, if someone could just point me to a copy of the sample, I would be totally fine with using that tool, too.


texts: Collapsing layer in current data frame of ArcMap using ArcPy

texts:
Online ArcGIS resource help provides a list of layer properties which can be used as needed. I am not able to find a way to collapse a layer using ArcPy in it.

Is a way to collapse a layer in the current data frame of ArcMap using ArcPy?