texts: Using Fiona to write a new shapefile from scratch

texts:
I am new to Fiona, and trying to wrap my head around reading a csv file and writing out an ESRI shapefile.  I have most of it down, but am missing two things; (1) defining the projection/spatial reference and (2) the attributes actually being populated with the values from the csv file.

So my questions are:
1. How do I define the spatial reference, either by defining an EPSG or a custom projection
2. Why is my data not writing as attributes?

Thanks in advance for any guidance.

import csv
from shapely.geometry import Point, mapping
from fiona import collection
from fiona.crs import from_epsg

schema = { 'geometry': 'Point',
 'properties': {'dx': 'float:13.3',
            'dy': 'float:13.3',
            'dline': 'str',
            'dtrace': 'int',
            'ddepth': 'float:9.4',
            'dtrash': 'float:9.4',
            'sx': 'float:13.3',
            'sy': 'float:13.3',
            'sline': 'str',
            'strace': 'int',
            'sdepth': 'float:9.4',
            'strash': 'float:9.4',
            'bx': 'float:13.3',
            'by': 'float:13.3',
            'bline': 'str',
            'btrace': 'int',
            'bdepth': 'float:9.4',
            'btrash': 'float:9.4',}}

with collection(
shapeout, "w", "ESRI Shapefile", schema) as output:
with open(finalout, 'rb') as f:
    reader = csv.DictReader(f)
    for row in reader:
        point = Point(float(row['dx']), float(row['dy']))
        output.write({
            'properties': {'dx': 'float:13.3',
            'dy': 'float:13.3',
            'dline': 'str',
            'dtrace': 'int',
            'ddepth': 'float:9.4',
            'dtrash': 'float:9.4',
            'sx': 'float:13.3',
            'sy': 'float:13.3',
            'sline': 'str',
            'strace': 'int',
            'sdepth': 'float:9.4',
            'strash': 'float:9.4',
            'bx': 'float:13.3',
            'by': 'float:13.3',
            'bline': 'str',
            'btrace': 'int',
            'bdepth': 'float:9.4',
            'btrash': 'float:9.4',},
            'geometry': mapping(point)
        })

With the code updated to the "newer" (aka, not outdated) process:

with fiona.open(shapeout, 'w', crs=from_epsg(3996), driver='ESRI Shapefile',
            schema=myschema) as output:
with open(finalout, 'rb') as f:
    reader = csv.DictReader(f)
    for row in reader:
        #geometry
        point = Point(float(row['dx']), float(row['dy']))
        #attributes
        prop = {'dx': float(row['dx']),
                'dy': float(row['dy']),
                'dline': str(row['dline']),
                'dtrace': float(row['dtrace']),
                'ddepth': float(row['ddepth']),
                'dtrash': float(row['dtrash'])},
        print{'geometry': mapping(point), 'properties':prop}

Gives me this:

{'geometry': {'type': 'Point', 'coordinates': (-245215.3333, -975226.3333)}, 'properties': ({'dtrace': 24387.0, 'ddepth': 1433.2458, 'dtrash': 1433.25, 'dx': -245215.3333, 'dy': -975226.3333, 'dline': '      lsl0904.201302'},)}
{'geometry': {'type': 'Point', 'coordinates': (-245215.0, -975251.0)}, 'properties': ({'dtrace': 24386.0, 'ddepth': 1432.7268, 'dtrash': 1432.73, 'dx': -245215.0, 'dy': -975251.0, 'dline': '      lsl0904.201302'},)}

But as soon as I change to:

output.write({'geometry': mapping(point), 'properties':prop})

I get the error message:
AttributeError: 'tuple' object has no attribute 'keys'


texts: Configure styling of OSM with T-Rex vector tile server

texts:
I am using T-Rex as vector tileserver, osm2pgsql for importing the .osm.pbf file into the PostGIS database (table format - planet_osm_line, planet_osm, point, planet_osm_polygon, planet_osm_road).

I want to display the osm vector tiles with style of osm-bright.

T-Rex is using .toml configuration file, where are specified layers with sql queries.

Do you have such a pre-configured style configuration for mentioned table structure ?


texts:
I have installed QGIS 2.4.0 on Windows. I also downloaded an OSM file form Geofrabrik (SouthAmerica).

I have been searching an option (plugin also) in QGIS for importing the OSM file to my PostGIS with no results.

Is there an option in QGIS or will I have to do it another way? (maybe with osm2pgsql manually)


texts: Where can I obtain 3D data for Singapore?

texts:
I am looking for providers of 3D maps/models of Singapore.

This might be either from photogrammetry/laser techniques or based on a DEM + building footprints/landuse with a height attribute.

Is any such data available either in public repositories or from specialized companies?


texts: Excel to GDB field mapping and importing Excel data in geodatabase?

texts:
I've got one Excel file which consists of data. I need to export that data to a gdb table. There are challenges. The gdb schema does not match with the Excel schema (I mean column names are different in excel and gdb table). I have to map each columns of the GDB with Excel.

I do not need to write the code.

I need myself if I can do it manually using the toolbox tools.

I am very new in Toolbox handling and arcpy as well.


texts: Aggregating polygons with space between them

texts:
Im currently editing a large dataset that has large chunks removed throughout (buffered road network). The areas adjacent from each side have coinciding attributes, but have this large erased area in between them as shown below:

The blue dotted lines below is where i'm looking to fill in the polygons for example.

Then by some workflow I need the final area to look like the following:

Is there a potential workflow (or tool combination) I could use to almost "reverse erase" this?


texts: Clipping raster with polygon mask but keeping outer edge cells of mask in QGIS?

texts:
In QGIS, how to clip a raster with a polygon mask, but keep the outer edge cells of the raster touching the mask?

In other words: how to create a one raster cell outer buffer for the masking polygon?


texts: Extract latitudes and longitudes from a government database

texts:
I'm trying to get the latitudes and longitudes of various data points on the Alberta Water License viewer. The mapping tool has the points displayed on it, and I'm able to click on them and get more information, but the only way that I can get the Lat/Long of the points currently is to open Google Earth and try and match the location exactly and manually write out the Lat/Long from that.

Is there some way I can extract this information more easily?

I don't have much GIS experience.


texts: Error in converting values in the field calculator in QGIS

texts:
I would like to calculate the flux of water in a river system.
The flux is given in km3/year and I would like to convert it to m3/month. My syntax is:

"flux"*1000000000/12

If I execute this the following error is shown:

  An error occurred during the evaluation of the calculation expression:
  Can not convert ',' to floating-point number.

The same error I keep getting if the to_real function is running.

Does somebody have an idea?


texts:
I have experience working in python but limited knowledge of Java Script.
I need help loading in a external Geojson located on my computer into a leaflet map.  I read that I need to use a technique called JSONP in order to bypass cross-domain securities.  This is where I am stuck.  I got a map working by hard-coding a single data point, into the map, but I still have hundreds more so that is not a repeatable and viable option.
Data
This is the URL form of the geoJSON.  I don't explicitly want to use the URL form as this is just a sample data set as I work through this problem.  I have another data set that is only in JSON form.  I will have to convert it to geoJSON before I can use in leaflet.  If I am wrong in this knowledge, loading a json file into leaflet is my question.
geojson
http://oceanview.pfeg.noaa.gov/erddap/tabledap/erdCalCOFIcufes.geoJson?longitude%2Clatitude%2Csardine_eggs&cruise=%22201504%22&sardine_eggs%3E=0&.draw=markers&.marker=5%7C5&.color=0x000000&.colorBar=%7C%7C%7C%7C%7C&.bgColor=0xffccccff
optional json (bringing this into leaflet is my end goal if it is possible)
https://coastwatch.pfeg.noaa.gov/erddap/griddap/jplMURSST41mday.json?sst[(2017-02-16):1:(2017-02-16T00:00:00Z)][(32):1:(38)][(-126):1:(-116)]
This is my working code
It has the map pulled up with the single point loading in as described above.

  <title>A Leaflet map!</title>
  <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.7.3/leaflet.css"/>

I have tried several tutorials on jsonp, but I just don't quite get the process yet.


texts: Showing only ONE layer from group layer in Leaflet?

texts:
I have the following code to solve my workflow. But it does not work properly.

The problem of the map behaviour appears to be within the map.getZoom function at the end of the code. That function should only be working when train layer is toggled on. But for now the getZoom is reacting on zoom even the train layer is not activated.

For clarification: there are two train layers which actually are the same layers coming from the same resource both. The difference is that aach of them is linked to a different style. Depending on zoom level one layer is displayed and the other is not. The both are grouped to a LayerGroup to have them both displayed as one layer in the LayerControl.

I am not a Developer.

Desired performance:

Opening the map with basemap active
Panning/zooming around
Toggle on the train layer in LayerControl
Train layer icons are displayed differently depending on the map's
zoom level

Actual performance:

Opening the map with basemap active
Zooming leads to displaying the train layer automatically --> Note:
displays without layer being toggled on in LayerControl
If layer is toggled on via LayerControl it shows iconTrain1 +
iconTrain2 --> two icons visible
Toggling off the layer in LayerControl removes the icons, but
zooming in/out displays them again (again without being toggled on
in LayerControl)


texts: Copy the layer_styles table from one database to another?

texts:
I am trying to copy the layer_styles table to another database. A quick and painless process is preferred, to make the movement of the layer_styles tables easy to repeat in the future.

I have tried two methods so far:
1. pg_dump and psql
2. Creating a dblink

Neither seem to work. The culprit seems to be the xml type in the layer_styles table.

I've had success using pg-dump and restore with other tables, however not the layer_styles table. Using pg_restore returns the following error:

   ERROR:  invalid XML content
    DETAIL:  line 1: StartTag: invalid element name
    <!DOCTYPE qgis PUBLIC 'http://mrcc.com/qgis.dtd' 'SYSTEM'>
     ^
    CONTEXT:  COPY layer_styles, line 1, column styleqml: "<!DOCTYPE qgis PUBLIC 'http://mrcc.com/qgis.dtd' 'SYSTEM'>
    <qgis simplifyLocal="1" maxScale="0" minS..."
     setval

The end result is an empty layer_styles table in the target database.

After creating the dblink, I run the following code:

CREATE TABLE layer_styles
AS
SELECT
from dblink('host=localhost
              user=postgres
              password=postgres
              dbname=database3',
              'select *
              from layer_styles') AS linktable(
        a integer,
        b varchar,
        c varchar,
        d varchar,
        e text,
        f xml,
        g xml,
        h boolean,
        i text,
        j varchar,
        k xml,
        l timestamp);

This returns the following error:

>  ERROR:  invalid XML content
>     DETAIL:  line 1: StartTag: invalid element name
>     <!DOCTYPE qgis PUBLIC 'http://mrcc.com/qgis.dtd' 'SYSTEM'>
>      ^

I have had previous success with the dblink when using different tables without an XML type.

Does anyone have any advice?


texts: Two Selections from one Layer in Processing Modeller

texts: I have the following problem: I try to do two parallel selections. For this I have applied before the "extract by location" algorithm. Now I would like to execute in the next step from this layer two parallel selections with the "select by expression" algorithm. However, this probably does not work. The modeller created me the two selection layers, but contain both the same records. So the one selection is carried out correctly, but for the other selection it takes the same value.


texts: Data Reviewer - how to check for polygon must be covered by line?

texts:
I am a new Data Reviewer user and have spent too much time trying to solve this so I knew it was time to ask for help. I have a polygon layer and a line layer.

They are in different feature datasets, so topology won't work here. I need the Data Reviewer equivalent to polygon must be covered by line so that I can find where the polygon and line are not lined up properly.

I've been trying different options in Geometry on Geometry but nothing has worked.

Can anybody suggest something to try, or point me to an online resource, as I haven't found anything relevant.


texts: Merge shapefiles to one tool not working?

texts: When using the data management tool "Merge shapefiles to one" in QGIS why is it copying one of my shapefiles instead of merging the two together as it should?


texts: Can't modify the layer using a script tool on the current map

texts:
I want a tool that does two things: Update symbology, and hide some fields. I created a script tool for that:

import arcpy

in_features = arcpy.GetParameter(0)
ref_lyrx = r"C:\ArcGIS\LAYER\Symbologylayer.lyrx"

aprx = arcpy.mp.ArcGISProject("CURRENT")
myMap = aprx.listMaps('Map')[0]
myLayer = myMap.listLayers("inputLayer")[0]

desc = arcpy.Describe(myLayer)
field_info = desc.fieldInfo
fieldsToHide = "HideThisField"

for i in range(0, field_info.count):
    if field_info.getFieldName(i) in fieldsToHide:
        field_info.setVisible(i, "HIDDEN")
        myLayer.transparency = 50
        arcpy.AddMessage(field_info.getVisible(i))
        arcpy.AddMessage(field_info.getFieldName(i))
        arcpy.ApplySymbologyFromLayer_management(myLayer, ref_lyrx)

aprx.save()

This code runs without any issues. The final for loop there is where the action happens:

Set field to hidden - But it does not hide on the attribute table, it is still there
Change layer transparency - Works, layer gets transparent
Print the status of the field - Prints "HIDDEN". If I delete the setVisible above, it shows "VISIBLE" instead.
Print the field name - The correct field that should be hidden shows here. Just for diagnosis.
Apply Symbol from another Layer - Does not work

The results are confusing. Why the transparency changes, but the rest does not? Why it change the status of my field to hidden, but this change does not reflect on map? Why ApplySymbology don't work (if I run the arcgis tool, it does).


texts: Projecting Microstation .dgn that has no coordinate system

texts: I have a set of .dgn drawings that have been sent to me from a client.  All of which were never assigned to a coordinate system upon creating them.  I am attempting to assign a coordinate system to this .dgn data and export to AutoCAD which I will then export to ArcMap.  Is this process possible?  Is it possible to Export Microstation .dgn files straight to a shapefile?


texts:
I'm trying to make my legend horizontal em leaftlet by I'm not sure what I have to change. My base legend is based on this tutorial. Ultimately what I would like to have is something like this. Any ideas? Is it just a matter of change the css?

Best,


texts:
I'd like to create a (web)map, that is showing the range of an Airplane around some airports.

I tried to calculate a buffer with the range of the airplane.
Here you can see the result here.

But now I have realized that the result is wrong because planes do not take the straight route but fly a curve because its shorter.

Is there a way to calculate the range with the shorter curve?


texts: How to calculate multiple distances between set of points and vector grid in PostGIS?

texts:
My set contains the following input data: A set of points and a hexagonal (vector) grid.

For each cell in the grid, I want to calculate the driving-distance to each point. I thought this would be a nice approach for calculation in a pgRouting enabled PostGIS database. At the end, I'd like to have the information stored in the attribute table of the Grid!

Unfortunately I'm not that much into PostGIS so I'm not sure about how to implement the algorithm correctly. What I did so far:

I have a street network that I loaded into PostGIS. I used pgr to create a topology. It seems to work!

[Edit 1]

I managed to simplify this a little. I am now using the heaxagon's centroids as points. I calculated the nearest network node for each centroid and wrote the source ID into a new table containing all centroids AND the source ID of the closest node from the network.

I did the same for the blue points so I have two tables to work with. I now want to route from each centroid to each blue node. Since I have very little understanding of PostGIS and Loops I thought someone could help me with the final step. The loop is meant to SELECT a centroid (source id) and route to all the blue point's source id. Then it should write the cost into the respective columnn. So It should look somewhat like this:

HEAXAGON ID | HEX SOURCE ID | COST A | COST B | COST C | COST D | ... | COST J
1           | 2341          | 250    | 600    | 564    | 993    | ... | 1200
2           | 6521          | 260    | 756    | 604    | 870    | ... | 1190
...         | ...           | ...    | ...    | ...    | ...    | ... | ...
4331        | 63525         | 1500   | 1874   | 120    | 451    | ... | 145

I have 10 Points and 4331 Hexagons, so this would need like 43310 Dijsktra Operations. I hope you can help me with an algorithm that solves this since I'm completely desperate right now!


texts: Using GeoXml3 to parse KML Markers loaded in Fusion Table Layer

texts:
I'm trying to apply the same code as this exmaple to zoom the google map to a KML Markers loaded into Fusion Table Layer.

The example code works for KML Polygon loaded into Fusion Table Layer but not for markers. Does anyone know how to apply that to KML markers?

kml example:

<Point><coordinates>-59.028645,50.772916,0.0</coordinates></Point>

script:

query = "SELECT * FROM  1eiOWq49sR8DXdymYALNdItE0uf-f4mNRbyJOVOU";
query.send(zoom2query);

function zoom2query(query) {
   // zoom and center map on query results
   //set the query using the parameter
   var queryText = encodeURIComponent(query);
   var query = new google.visualization.Query('http://www.google.com/fusiontables/gvizdata?tq='  + queryText);

    //set the callback function
    query.send(zoomTo);

      }

  function zoomTo(response) {
    if (!response) {
      alert('no response');
      return;
    }
    if (response.isError()) {
      alert('Error in query: ' + response.getMessage() + ' ' ++ response.getDetailedMessage());
      return;
    }

      FTresponse = response;
      //for more information on the response object, see the documentation
      //http://code.google.com/apis/visualization/documentation/reference.html#QueryResponse
      numRows = response.getDataTable().getNumberOfRows();
      numCols = response.getDataTable().getNumberOfColumns();
      var geoXml = new geoXML3.parser();
      var bounds = new google.maps.LatLngBounds();
      for (var i = 0; i < numCols; i++){
        if (FTresponse.getDataTable().getColumnLabel(i) == 'geometry') {
         var ColIndex = i;
        }
      }
      if (!ColIndex){
        alert('Geometry column "geometry" not found.')
      }
      for (var i = 0; i < numRows; i++){
        var bounds = new google.maps.LatLngBounds();
        for(i = 0; i < numRows; i++) {
            var kml = FTresponse.getDataTable().getValue(i,ColIndex);
            geoXml.parseKmlString("<Placemark>"+kml+"</Placemark>");
            bounds.union(geoXml.docs[i].bounds);
        }
        // zoom to the bounds
        map.fitBounds(bounds);
      }


texts: I am using Ubuntu 14.04, Postgresql 9.6. Previously I had working version of Postgis 2.3 but after upgrading to 2.4 when I query something from database with geometry objects I get ERROR:  could not access file "$libdir/postgis-2.4": No such file or directory. If I make count(*) on the same table, it works well. How to solve this problem?


texts: GeoTools example for Shapefile export to already existing PostGIS table

texts: I'm a newbie to GeoTool API, and I'm looking for a way to export the specific shapefile to the specific already existing PostGIS table with GeoTools (the underlying problem I face is that PostGIS tables are configured such as I can only view them with uDig but not other GIS; uDig can't edit them; I draw some additions on the separate layer, export then to the shape file and want to add them to the existing PostGIS tables whose data is used as a background when I draw the supplement in separate layer. Since uDig fails to export the shapefiles to the tables I desire, I have to do it manually, but I'd like to have some automation there)/


texts: Publish CSV to Enterprise Portal as table

texts:
I'm following this help doc.
https://geosaurus.maps.arcgis.com/home/item.html?id=a1db6db172bc49a8932daacc2ed3d3ac#preview
I'm able to add the CSV to the Portal, but I cannot get it to publish. I need it published as a table. It doesn't have coordinate fields and it's not intended to be a spatial feature.
from IPython.display import display
from arcgis.gis import GIS
import os
gis = GIS('Home')

csv_file = r'C:\Users\jpilbeam\Downloads\c19_Vaccine_Current.csv'
csv_item = gis.content.add({}, csv_file)#add csv to Portal

csv_layer = csv_item.publish(None, file_type='csv')

Error:
--------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
In  [15]:
Line 1:     csv_layer = csv_item.publish(None, file_type='csv')

File C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\lib\site-packages\arcgis\gis\__init__.py, in publish:
Line 10353: elif not buildInitialCache and ret[0]['type'].lower() == 'image service':

KeyError: 'type'
---------------------------------------------------------------------------

Other attempts:
csv_layer = csv_item.publish()
KeyError: 'type'

csv_layer = csv_item.publish('csv')
ValueError: dictionary update sequence element #0 has length 1; 2 is required

csv_layer = csv_item.publish(publishParameters={"type":"csv"})
TypeError: publish() got an unexpected keyword argument 'publishParameters'

I also tried what the help doc does. It has a text field with place name values specified in the address_fields parameter. My CountyName field is also a text field with place names. But, this also errors.
csv_flayer = csv_item.publish(None, {"Address":"CountyName"})
KeyError: 'type'


texts: Openlayers in QGIS - Errors in opening google satellite/

texts:
End goal: To put contour lines with lat&long layers from a google earth location (Rock Piles beach, HI) in QGIS. I plan to use CRS: EPSG:3760

Does the Plug in (Openlayers / QuickMapServices) already have the D.E.M (xyz file) or do I need to get it from somewhere else?

I installed both 'openlayers plugin' and 'quickmapservers' plugins in QGIS. Then when I go to the web tab > openlayer plugin > Google Maps > Google physical.
This message shows:

The google terrain still downloads, but when I double-click the layer to see the properties, this message pops out:

Is there a way to fix this?


texts:
How to solve the problem I show in both images. I have a centered map where I want, when I click to print I would like it to print almost from the screen and reproduce the same position in the print.

However, when I print the map it leaves in a completely different position from the one I want and with more zoom than normal.

What is the best solution to solve this problem? Attached is the two images of the map, one of how it is and another of how it is in the print. I wish they were both alike.

Resumming: How can centralize my map in the print as without print? Attached is the two images:  and . I need a tutorial or some how fix this, please. I hope you understood, sorry for my english.


texts: Using reducer over space and time in Google Earth Engine?

texts:
I am trying to use in Google Earth Engine a reducer (histogram) on an ImageCollection, computing it for multiple shapes for each image.

For a given image, I use reduceRegions(..., ee.Reducer.mean()) and this work, returning a feature collection with property mean.

However, when I try to apply this to all the image of the ImageCollection, using map(), it still returns a FeatureCollection, but the mean property is empty.

How can I solve this? Why is my result good on an image, but weird on an ImageCollection

var rec1 = ee.Geometry.Polygon([[[-89.78276, 40.160312], [-89.78293, 40.153424],
      [-89.77108, 40.153424], [-89.77117, 40.160377]]]);

var rec2 = ee.Geometry.Polygon([[[-89.782848, 40.15335], [-89.782848, 40.14607],
      [-89.774522, 40.14607],[-89.774265, 40.15322]]]);

var featCol = ee.FeatureCollection([ee.Feature(rec1), ee.Feature(rec2)]);

/// get ImageCollection, and Image
var CDL= ee.ImageCollection("USDA/NASS/CDL")
var CDL_2015 = ee.Image('USDA/NASS/CDL/2015');

/// For one image: works
var fieldStats_oneIm = CDL_2015.select("cropland").reduceRegions({
collection: featCol,
reducer: ee.Reducer.mean(),
scale: 30
});

/// For all images: does not work?
var fieldStats_ImCol = CDL.map(function(image) {return
    image.select("cropland").reduceRegions({
    collection: featCol,
    reducer: ee.Reducer.mean(),
    scale: 30
})});

print(fieldStats_oneIm, "fieldStats_oneIm")
print(fieldStats_ImCol, "fieldStats_ImCol")


texts: How can I import a GeoPDF file as a layer in QGIS 3?

texts: I want to open a USGS quadrangle (such as this) in QGIS 3, but I can't find documentation for how to do this.


texts: How to add N+1 second using ArcPy cursor for Time Field value?

texts:
I have a field that I would like to add one second on for each feature using field calculator and Python.  My field name is ISSUEDTIME_NEW I want to start my time at 00:00:00 and add one second on to each feature.

Here is what I have so far for a search cursor, however, this does not loop as expected and I want to order by another field, ISSUEDDATE so that my earliest time in my order field is my earliest hh:mm:ss in my ISSUEDTIME_NEW field.

i.e. since 9/10/2017 12:31:00 PM is my earliest time in my order field then it equals 00:00:00 in my ISSUEDTIME_NEW

What change should I make in the code to add 1 second onto 00:00:00 for ISSUEDTIMENEW and ordering by ISSUEDDATE?

import arcpy
import datetime

fc = 'fc'
sql = ('ORDER BY ISSUEDDATE')
field = ('ISSUEDTIME_NEW')

with arcpy.da.UpdateCursor(fc, field, sql) as cursor:
    for row in cursor:
        a = datetime.datetime(100, 1, 1, 12, 00, 00)
        b = a + datetime.timedelta(0,1) # days, seconds, then other fields.
        while row[0] is True:
            row [0] = b.time
        cursor.updateRow(row)`

RuntimeError: An invalid SQL statement was used. [SELECT ISSUEDTIME_NEW, OBJECTID FROM FC WHERE ORDER BY ISSUEDDATE]


texts: How can i use ihookhelper interface in python?

texts:
I use ArcObject in python to develope a tool. I import comtypes in python to access  the arcmap layers in python and use ifeaturelayer interface. I have to set featurelayer to phookhelper.focusmap.get_layer(i) to access the layers in ArcMap. The problem is that I can't find get_layer method in python. Is there any way to access the ArcMap layers in python using ArcObjects?

for (int i = 0; i < phookhelper.FocusMap.LayerCount; i++)
            {
                if (phookhelper.FocusMap.get_Layer(i).Name == comboBox1.Text)
                {
                    pfeaturelayer = (IFeatureLayer)phookhelper.FocusMap.get_Layer(i);

                    pfc = pfeaturelayer.FeatureClass;


texts: Add separator and text after numeric scalebar in QGIS 3

texts:
In the numeric scalebar in QGIS 3 is it possible to add a thousands seperator of ',' and also some text after the value?

This should be "1:500,000 at B0" for example.

I have tried Display and Font options but can't find anything to format the value as numerical with a thousands seperator or add a suffix text.


texts:
OGC filters not working on WFS service It only seems to serve up unfiltered data or a layer with a attributes table but no features as if nothing met the criteria. I have tried repeatedly with criteria that I know is in the dataset. Am I doing something wrong here?
Using Qgis 2.14.11 built in WFS client

http://openmaps.gov.bc.ca/geo/ows?version=1.3.0&SERVICE=WFS&VERSION=1.0.0&REQUEST=GetFeature&TYPENAME=pub:WH SE_FOREST_VEGETATION.VEG_COMP_LYR_R1_POLY&SRSNAME=EPSG:3005&FILTER= <ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
 <ogc:PropertyIsEqualTo>
  <ogc:PropertyName>MAP_ID</ogc:PropertyName>
  <ogc:Literal>082M099</ogc:Literal>
 </ogc:PropertyIsEqualTo>
</ogc:Filter>

When I run it with the filter removed from the URL, back to the CRS (EPSG:3005) it works but I usually run into a download limit long before getting the features I am after. If I run it with different criteria it works at it for what seems to be a appropriate amount of time considering the amount of data that It "should" be returning but no features. having the same problem with different datasets served up here.  I am suspecting that either the server just will not serve filtered data or the method I need to use is completely different somehow.


texts: Calculate field does not update using ArcPy

texts:
I want assign 33 if the field no_apt is zero and 34 if no_apt is any else value. The model run without any error, but no value are assigned the field GGV_kW. What is wrong with the code since it doesn't work?

The code used is:

in_Table=nodesFeatureoutall
fieldName="GGV_kW"
expression="GGV_calculation( !no_apt!, !GGV_kW!)"

codeblock = """

def GGV_calculation(no_apt,GGV_kW):

    if no_apt==0:
        return 33
    else:
        return 34"""

arcpy.CalculateField_management(in_Table,fieldName,expression,"PYTHON",codeblock)


texts: Finding whether a particular band value is consistently less over a period of time in an image collection

texts:
The code below gives me forest regions of India from the Hansen data set and then I have filtered the landsat images of these regions from 2013-2019 on a quarterly basis. My output image is such that each pixel has its value as the number of pixels with NDVI less than 0.25 at that particular pixel location over these years.

Let's suppose, in the output image, one of the pixel has a value = 6, that is, at that particular pixel location, there are 6 pixels over these years which had NDVI < 0.25. I have to find out whether these 6 pixels occurred consecutively or not. That is, I want to know whether NDVI was constantly low in the consecutive quarters and years. If that is the case, then I can conclude that region is under threat. I also want to know the time period when the NDVI was consistently low.
Is there a way I can do this?

// FILTERING OUT FOREST COVER IN INDIA

var treeCanopyCover = hansen.clip(india).select('treecover2000');
var greater25 = treeCanopyCover.gte(25);

// NDVI FUNCTION

function addNDVI(image) {
    var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');
    return image.addBands(ndvi);
    }

// FILTERING THE COLLECTION

var startDate = ee.Date('2013-01-01');
var endDate = ee.Date('2019-06-30');

var filtered = l8.filterBounds(indiaforest2)
              .filterDate(startDate, endDate)
              .filterMetadata('CLOUD_COVER','less_than',25)
              .map(addNDVI);

//print(filtered);

var maskedCollection = filtered.map(function(image){
  var x = image.updateMask(greater25);   // filtering over forest regions
  return x;
});

// filtering quarterly wise

var map_m = function(y)
{
  var months = [1, 5, 9];
  var filtered_col = months.map(function(m)
  {
    var start = ee.Date.fromYMD(y, m, 1)
    var end = start.advance(3, 'month');
    var filt = maskedCollection.filterDate(start,end).mean();

    return filt.set('year', y)
              .set('month', m)
              .set('system:time_start',ee.Date.fromYMD(y,m,1)) ;
  });

  return filtered_col;
};

var years = ee.List.sequence(2013, 2019);
var quarterlyimages = years.map(map_m).flatten();
print(quarterlyimages);

var z = ee.ImageCollection(quarterlyimages);

// discarding null images

var nullimages = z.map(function(image)
{
      return image.set('count', image.bandNames().length());
})
    .filter(ee.Filter.eq('count', 13));

// counting no. of pixels with NDVI < 0.25

var NDVI_count = nullimages.select('NDVI')
  .map(function(img) { return img.updateMask(img.select('NDVI').lt(0.25)) })
  .count();

print(NDVI_count);

var ndviViz = {min: 1, max: 20, palette: ['green', 'yellow', 'orange', 'red','violet','indigo', 'blue']};

Map.addLayer(NDVI_count, ndviViz);


texts: How to find the intersection point of polygon geometry and multilinestring geometry in postgis

texts: I have a polygon geometry and multilinestring geometry. How do i find the intersection of both these geometries. Both the geometries are in same reference system. Does st_intersection gives me the intersection point???


texts:
I have a dxf file which was converted from a CAD file. I used QGIS to convert it to a shp file...
This is when I realized that the text associated with the drawing is actually a separate layer...
What I need is one layer where the text is in the attribute table.

What would best the best way to accomplish this?


texts:
New to Python.

I have an mxd file file with three layers and one is a layer of elementary school boundaries.

Goal: Make a pdf doc with a page for each elementary boundary.

arcpy.SelectLayerByAttribute_management(Counties, "NEW_SELECTION", "Boundary = 'PineCone'")
ExportToPDF
Next boundary polygon
Append PDF
Save

My question: Do I need to make a physical file of the elementary boundary layer before I can “SelectLayerByAttribute” for each boundary? Or is there a way to loop through the boundaries by making the selections straight from the mxd file.Below is not complete code, just the thought process.

import arcpy
from arcpy.mapping import *
mapDoc = MapDocument(r"C:temp\boundary.mxd")
pdfElementary = r”C:\Temp\Elem_Bndries”

???
???

for bname in boundaryList
arcpy.SelectLayerByAttribute_management(Counties, "NEW_SELECTION", "boundary = bname")
ExportToPDF (mapDoc, pdfElementary, "PAGE_LAYOUT")
pdfFinalBoundary  = appendPages(pdfElementary)

pdfFinalBoundary .saveAndClose ()
del pdfFinalBoundary


texts: Elevation data from 3D shapefile of buildings to raster DEM

texts:
I have acquired LoD2 building models in the 3DShape format (Polygon ZM). This LoD2 data needs to be prepared for a R-script of solar irradiation calculations on roof surfaces, but which can only process raster data.

Is there a way to 'rasterize' the 3D Shapefile into a digital surface model (DSM)?

Lets say for example, that I want a spatial resolution of 0.25 x 0.25 meters for each raster cell, than the height of the highest feature (i.e. the roof polygon) should be extracted from the LoD2 data at the location where it intersects with the center of each raster cell. So that in the end I get a 'DSM' of the shapefile.

Sample LoD2 data can be downloaded here: http://vermessung.bayern.de/file/zip/6817/LoD2_3DShape.zip

I have found nothing in the ArcGIS help resource center, and no GRASS function either. ArcGIS, QGIS, GRASS and R are available to me.

If you know a alternative way to calculate solar irradiation on a 3DShape surface, you would be of help as well. I have irradiation measurements for an fixed normal (=horizontal) surface in the area. What my R script does is it calculates the effect of the rooftop geometry and the shadows of nearby objects (= other buildings) on the measured irradiation at 1 hour time-steps over an entire year and then summarizes the W/m² to evaluate which rooftops get the most sunlight.


texts: Exporting Vector Layer with OGR2Layers Plugin?

texts:
I'm using QGIS 2.0 Dufour on my Linux (Ubuntu 12.10).I'm download shapefile
of Croatia country from this site (http://www.gadm.org/country) and i added a new vector layer where i created some polygons. (a picture is attached).
When i created OpenLayers HTML map using OGR2 layer plugin, desired polygons are not displayed in open street map.

Where am I wrong, if it is even possible?

Vector layer in QGIS


texts: Some isolated point symbols not rotating while rest do using QGIS

texts:
I'm trying to rotate some svg geology structure symbols in QGIS, using point data originally from a geoJSON file which I converted to geopackage.
I used the expression "attitude_strike"-90 to apply the proper rotation to the symbols, and it almost works great. About 95% of the symbols rotate as expected, but a small number of them do not rotate.
These are symbols of one data type, with the same symbol, etc.
Any ideas what could be causing this?
Here are some screenshots, first showing the various unrotated symbols, then showing the attribute table for one of those symbols, which has an strike of 10.


texts: ArcGIS Explorer iOS App can't load map with bing basemap

texts: In the Explorer iOS app, I am getting an error that says "Missing bing maps key" when opening a map with a bing map as the basemap. The bing basemap works fine on the full website and mobile website and I have verified that my bing maps key is still valid. This has been happening for a few weeks now. What is going on?


texts: QGIS and my Network Attached Storage Server don't play well together

texts:
I use QGIS on a desktop and multiple laptops. I store my project files on all the local computers. The project files are synced to a Synology NAS server using "Synology Drive Client" software.

I have been using this setup for years using ArcGIS with no issues, however since using QGIS, sometimes my updated edited geopackage files get overwritten by an old version and I have lost all my edits. My colleague working on the same project is having the same issue.

Does anyone have any solutions?


texts: How can I select all streams downstream of a polygon

texts:
I have an ordered stream network and I want to select all stream features that contain water that has come from a certain area. To start with I'm running a select by location query to find all stream segments that cross my polygon. The problem is that this does not select the stream segments further down stream that these segments flow into.

What I need to do is create a selection that will select stream segments that touch existing selected segments only if they have a higher number for the stream order attribute, then repeat this until the result does not change. I expect I have to do this in python, but I'm not sure of a few things: how to use the select by location and  select by attribute at the same time (select intersecting lines that have a value above a threshold), how to make the threshold relate to the line you're checking the intersection with (select intersecting lines that have a higher value than the line they're intersecting with, how to get this to iterate through until the selection isn't changing.

This feel quite complicated, but I'm sure people will have run into this problem before so there might even be an easier way, or a tool. Is there?

This image is a manual section that shows what I'm trying to achieve - I want to  select the streams that start in the green area and everything they flow into.


texts: MMQGIS: delimiter error on CSV import for merging

texts:
I have a shape of my city and I have a CSV file, containing attributes I'd like to join into the shape of the city. In selection of the CSV file, I get the following error message by MMQGIS.

Source of the CSV: public server of election data. Things I tried to fing the error:

converted the file to UTF-8 encoding
converted it to ASCII
replaced \n with \r\n
replaced ; delimiters with | and ,
replaces decimal characters with . and ,
open and resave in Libre Office
open and resave in Notepad++
use CSV lint for verification
view the file in Libre Office to compare ammount of columns
split the CSV (~30 lines) to merge them step by step to find broken lines
found the line where the import seems to break
still can't find whats broken in this line

Is it possible that this tool is causing the error itself?

Update: I just skipped the first 5 lines and tried the next 5 lines. Works without error. Maybe it is the length of the lines? Or the total ammount of data in the CSV? That would be odd because it is not much data.

If I merge only a few lines, all other lines with no new CSV attributes disappear from my shapefile.

Shapefile: https://www.offenedaten-wuppertal.de/dataset/40c4cc39-08d2-44df-9257-94edc7af3578/resource/40c4cc39-08d2-44df-9257-94edc7af3578

Example data to merge: (Source: https://wahlen.wuppertal.de/BW2017/05124000/html5/Bundestagswahl355.csv)

  06;friedrichsberg;7;7;7150;4912;68.7;4834;78;4855;57;1321;27.33;1182;24.35;1635;33.82;1337;27.54;349;7.22;406;8.36;475;9.83;598;12.32;349;7.22;559;11.51;526;10.88;520;10.71;59;1.22;37;0.76;0;0;15;0.31;114;2.36;58;1.19;0;0;12;0.25;0;0;3;0.06;0;0;8;0.16;6;0.12;4;0.08;0;0;1;0.02;0;0;30;0.62;0;0;8;0.16;0;0;9;0.19;0;0;4;0.08;0;0;10;0.21;0;0;4;0.08;0;0;1;0.02;0;0;44;0.91;0;0;5;0.1;
  11;brill-arrenberg;6;6;6141;4697;76.49;4659;38;4670;27;1319;28.31;1155;24.73;1433;30.76;1001;21.43;542;11.63;638;13.66;473;10.15;621;13.3;435;9.34;704;15.07;313;6.72;326;6.98;31;0.67;22;0.47;0;0;4;0.09;105;2.25;70;1.5;0;0;5;0.11;0;0;3;0.06;0;0;4;0.09;8;0.17;7;0.15;0;0;0;0;0;0;21;0.45;0;0;11;0.24;0;0;17;0.36;0;0;3;0.06;0;0;20;0.43;0;0;5;0.11;0;0;4;0.09;0;0;27;0.58;0;0;2;0.04;


texts: How to extract NDVI values into polygon vector map with QGIS

texts: I have calculated the NDVI values of a satellite image using the raster calculator. Now I need to have these NDVI values extracted into the attribute table of my polygon vector file. Is there a plug-in that would help? Thanks!


texts: I have a need to to convert a CSV file that was surveyed in NAD83 Zone 12 to NAD83 Zone 13 coordinates, is there an online tool that can do this for me?


texts: Note QGIS raster calculator defaults to not being accurate

texts:
[Edited title of post based on fact that @wingnut solved it]
I'm trying to normalise a population density map such as this one (e.g. to a maximum of 255). The maximum value in the file to begin with is 27500.1. However, if I enter ("gbr_pd_2020_1km@1" * 255) / 27500.1 as a formula in the Raster Calculator it only gives me values up to 198.
What makes me further think it's a bug is if I just enter as a formula "gbr_pd_2020_1km@1" / 27500.1, it should give me values up to 1, but it only gives me values up to 0.8. And even if I just enter itself * 1 it only gives me values up to 21367!!! Please see

What I think must be happening is that the 27500.1 value is only in one pixel (which represents one area in London) and somehow the raster calculator is not looking at every pixel. But I can't see any settings to change this.
I've also tried doing the same via 'Processing Toolbox' then GDAL Raster Calculator, but it gives the same incorrect results.
I've tried this in QGIS 3.16.5 and 3.18 (the latest LTR and beta releases available), and am running it on the latest version of Mac OS (Big Sur 11.2.3).


texts: QGIS 3: How to filter OpenStreetMap nodes by their role in a relation

texts:
I want to export a relation from OSM, style it and export it again to SVG or some other format to use it in other applications.

For fetching the data I use QuickOSM and the following query:

<osm-script output="xml" output-config="" timeout="25">
  <union into="_">
    <id-query type="relation" ref="9022117"/>
    <recurse from="relations" type="relation-node"/>
    <recurse from="relations" type="relation-way"/>
  </union>
  <union into="_">
    <item from="_" into="_"/>
    <recurse from="_" into="_" type="down"/>
  </union>
  <print e="" from="_" geometry="skeleton" ids="yes" limit="" mode="body" n="" order="id" s="" w=""/>
</osm-script>

In my example I am using relation 9022117, a train connection. I want to map the connection as a red line with the stops. I got the data, the red line and I got the nodes styled as I want. However right now all of the nodes are styled that way.

When looking at the relation I can see the member nodes like this:

...
    <member type="node" ref="6076636356" role="stop"/>
    <member type="node" ref="6076650267" role="stop"/>
    <member type="node" ref="6077173089" role="stop"/>
    <member type="node" ref="6076666707" role="stop"/>
...

How can I filter the nodes in QGIS, so that I can only style/display the nodes which are set as members with role stop in the relation? QuickOSM seems to not include the members.


texts: Using gdal_translate for multiband raster to convert into .XYZ file?

texts:
I have raster containing 2 bands (NDVI and DEM). This raster was created in SNAP using "Add Elevation Band" option, "DEM Model:SRTM 1sec Grid" to the NDVI product I had generated earlier.
I could not find any option in SNAP to make the conversion of this combined product into an ASCII file. There is however one option in SNAP to export the product to a CSV but the result of that process were not as expected. So I switched to QGIS. Here I tried the Raster>Conversion>Translate option, which ultimately uses the gdal_translate tool.

So in QGIS-Translate I could convert the .tif to .XYZ(ascii) but the output had only one band value along with the lat-long for each pixel. I got the warning "Warning 6: XYZ driver only uses the first band of the dataset."

I used the option -b 1 -b 3 to include the bands I want in my output file but still the same warning kept coming.

My gdal_query looks as below,

gdal_translate -of XYZ test-op-1p-ndvi_DEM.tif gdaltranslate-op1DEM.xyz -b 1 -b 3

I tried using the other option in this query like -ot with Byte format (as mentioned in one related question here in GIS-StackExchange) and also by removing the -of XYZ parameter

I have also tried gdal_translate through the OSGeo4W desktop tool but got same warning there too,

After searching hard enough, I could not find anything satisfactory. So now I am trying the long way by generating 2 separate .XYZ for each band and then combine through python script, parsing and combining the 2 columns.

Can anyone please provide any hint or suggestion. Even if someone can link to similar questions in this forum would be really helpful, as i have already tried solutions from many related questions but could have missed something important.

Screenshot of .XYZ with only NDVI values

Screenshot of .XYZ with only DEM values


texts:
I need to georeference this .tiff. I see coordinate values ranging South-North 550000-554000 and West-East 10000-12000. Do I need to know the coordinate system to georeference the image
?


texts:
I run QGIS 3.18.1 on Ubuntu 20.04, I miss the newsfeed on the welcome page:

Newsfeed is activated in the settings:

I was able to solve the missing Recent Projects according to the answer on this issue: QGIS 3.10 Recent Projects List missing but I was not able to access the newsfeed the same way. Any hints?


texts: QGIS 2.6 automatically zooms to layer extent

texts: Since I upgraded to QGIS 2.6 a mouse click (left or right) on any layer returns in a automatic zoom to the layers extent. Does anyone know how to switch off this function?


texts: Concurrent users using a DLL to access ArcObjects?

texts:
So currently I have some ASP.net code where users are updated GIS through a DLL file that uses ArcObjects to talk to the MapService and layer.

My question is, will there be concurrent user issues? What if 10 people are submitting updates to the same items that need updating in ArcGIS map service? Will the DLL freeze up? Will ArcGIS service freeze up?

Or does it handle concurrent users in a Queue-style via the ArcGIS service, the DLL, or IIS itself maybe realizes the DLL is in use?


texts:
I recently upgraded QGis on Win 8.1 to 2.18.10 and as I opened a project file I had created before upgrading, there seems to be things that do not work.

Can't select any menu item. Hovering over them with the mouse doesn't highlight them and when I click on them nothing happens
when I try to select layers from the layers side panel, nothing happens
when I click on toolbar buttons, the drawing zooms out....
even if I do not open an old project, and try to create a new one QGis is pretty unusable like this

I see that there is a newer release. I wonder if it addresses problems of this sort or can it be that the installation/upgrading process failed somehow?


texts: Showing Statistics for Multiple Attribute using 'Show Statistical Summary' tool in QGIS

texts:
I need on the fly statistic as I select point on canvas. What the tool does though is give you the statistic and the given value for the Attribute selected.

What if I need the value of multiples Attributes. let's say I want the sum of an attribute named Area, and of an other attribute named perimeter. I tried using the equation tool to concatenate the data like so:

  'area:' || area || 'perimeter:' || perimeter

But that doesn't create a new column for new different attributes and now the value just comes back empty.


texts: Attachments not saving in an Enterprise Geodatabase

texts: I have an ArcGIS Server 10.2.2 feature service.  The feature classes are saved in a SQL Server 2012 Enterprise Geodatabase.  The feature classes have attachments enabled.  They are non-versioned and have archiving enabled to allow syncing and offline editing.  I've added the feature service to ArcGIS Online.  I can edit the feature classes and add attachments with no errors.  But when I close the editor window and then identify the feature it says 'No Attachments Found'.  When I look at the geodatabase on the server there is a record in the __ATTACH table but when I Identify a feature it doesn't show any attachments.  Has anyone successfully setup attachments on a 10.2.2 feature service in an Enterprise Geodatabase?


texts: Calculate driving distance and time from 2 lat/lon points

texts:
I'm researching a way to calculate driving time and distance for 2 latitude and longitude points?

I'm unable to use any 3rd party controls or maps.  If possible it has to all be generated locally to my machine.  I'm working with very restricted data.

Can ArcGIS do this and if so what do I need to do this in a batch?


texts: Merging rasters in R - errors using mosaic function

texts:
I'm trying to use the mosaic function (package: raster) in R to join two geo-tiff rasters together. Both are derived from LandSat tiles. I understand that mosaic requires

  All objects must have the same origin, resolution, and coordinate
  reference system.

Unfortunately, the two tiles have different CRS

+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

and

+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

leading to

Error in compareRaster(x, extent = FALSE, rowcol = FALSE, orig = TRUE,  :
  different CRS

So, I converted them to the same crs by using projectRaster - I've tried reprojecting both rasters to both epsg:3857 and epsg:4326

However, having done this, in both cases when I run the mosaic function I get the error

Error in compareRaster(x, extent = FALSE, rowcol = FALSE, orig = TRUE,  :
  different origin

How can I fix this? Is there a different crs I could use that would use the same origin for both?


texts: What alternatives exist for EOBrowser.exe & LibraryLocator.exe for ArcGIS Desktop 10?

texts:
Original Question: Where are the EOBrowser.exe & LibraryLocator.exe tools in the ArcGIS Desktop 10 > Developer Kit?

I went to use LibraryLocator.exe and couldn't find it now that I'm on ArcGIS Desktop 10.  And that's when I realized there is no ArcGIS\DeveloperKit\Tools\ installed.

Do these exist in 10.0? Are they downloadable?  Or is expected that we use the website (which is what I generally do anyway)?  The fate of my SlickRun magic word for LibraryLocator.exe will be decided by the answers to this question.


texts: Spatial Join between lines and points does not work in QGIS

texts:
Here is the workflow:

I've imported lines from a shapefile
I've used the Random points along line tool to generate points all across the previously imported lines
Then in Vector > Data Management > Join Attributes by Location, I select both lines and points layers, select "Use first localized entity" and run the job

In the end, it never generates any points, as if it was unable to spatially find on which lines a given point is located (which is surprising since the points were generated from lines).

Am I missing something?
Tested in the latest version of QGIS and in version 2.18.

Is there any way to spatially join attributes of two layers together?

I've tested in ArcGIS Desktop 10 and it works fine.


texts: Nearest neighbour joins for two geodataframes using KDtree from scipy/sklearn

texts:
I have a number of large geodataframes and want to automate the implementation of a Nearest Neighbour function using a KDtree for more efficient processing. The process I want to achieve here is to find the nearest neighbour to a point in one dataframe (gdA) and attach a single attribute value from this nearest neighbour in gdB.

This is effectively the function performed by the NNjoin package in QGIS. The code written below (with help from this link GeoPandas: Find nearest point in other dataframe) succeeds in generating the same distances as NNjoin, but unfortunately does not return anything like the same attribute values.

OPTION 1

from scipy.spatial import cKDTree

def ckdnearest(gdA, gdB, bcol):
nA = np.array(list(zip(gdA.geometry.x, gdA.geometry.y)) )
nB = np.array(list(zip(gdB.geometry.x, gdB.geometry.y)) )
btree = cKDTree(nB)
dist, idx = btree.query(nA,k=1)
df = pd.DataFrame.from_dict({'distance': dist.astype(int),
                             'bcol' : gdB.loc[idx, bcol].values })
return df

SW_closest = ckdnearest(gdA, gdB, 'attribute_value')

OPTION 2

from sklearn.neighbors import KDTree

def ckdnearest(gdA, gdB, bcol):
nA = np.array(list(zip(gdA.geometry.x, gdA.geometry.y)) )
nB = np.array(list(zip(gdB.geometry.x, gdB.geometry.y)) )
btree = KDTree(nB)
dist, idx = btree.query(nA, k=1)
idx = idx[:,0]
index = idx.flatten()
df = pd.DataFrame.from_dict({'distance': dist.astype(int),
                             'bcol' : gdB.loc[index, bcol].values})
return df

SW_closest2 = ckdnearest(gdA, gdB, 'attribute_value')}


texts: How can I summarize all the pixels from a single raster using Google Earth Engine?

texts:
I am trying to summarize all the values from one single raster to get a value (with the total). I ve used this code and I get the total value of the summarize from all the pixels in the region but I don't know how to use that value for other expressions. Could anyone help me?

var image = ee.Image('LANDSAT/LC8_L1T/LC80440342014077LGN00')
.select(['B4']);

var region = geometry1

var sumB4 = image.select('B4').reduceRegion({
 reducer: ee.Reducer.sum(),
 geometry: geometry1,
 scale: 30,
 maxPixels: 1e9
 });

The total value is 654789874.465978 and now I need to use this value to introduce in a NDVI formula.

How can I do it?


texts:
I am trying to virtually put a blanket over a DEM so to say, in more technical terms, a convex hull over 3d point data. To get some insight in what I'm doing exactly, I will first explain my entire plan so it becomes clear what the purpose of the convex hull is.

I am trying to detect sand bars close to the shore of the Netherlands. I have point measurements of the depth and may plan is as follows:
first I remove the slope of the coast, secondly I invert the DEM and lastly I use ArcGIS' hydrology toolset to find the spots in the landscape where water will accumulate, these areas are my sand bars.

Here is a sample of my point data. It's transects from a part from the coast. And here is what I would like to to visually:

The order of the three images is arbitrary. I realised after making this image, it's also possible to first remove the angle and than invert it.

My question right now is that in order to remove the slope, my plan is to create a convex hull over the DEM, which would look like putting a blanket over the 3d landscape (enclosed on the bottom, but this is not relevant for me). After this I will subtract the point measurements form this straight slope, resulting in a flat landscape with indentations in it. Does anybody know of a python, arcgis or qgis library or tool that could do this? I have been looking at ArcGIs' minimum bounding geometry, but this only works for 2d data. R's Convhulln might work, but there might work, but I'd rather find something in Python.


texts: PyQGIS refresh layer tree programmatically

texts:
When programmatically applying a subset string to a vector layer (here: ax_flurstueck with object count 34265)

via

l = iface.activeLayer()
l.setSubsetString('gemarkungsnummer in ('2760','2765','2771','2780','2786','2800','2811','2817','2836')')

the map canvas refreshes automatically, but not the object count in the layer tree (displays still 34265 objects):

When moving the main window around or changing the width of the layer tree or something like that, the layer tree refreshes (21581 objects).

The question is, how to achieve this programmatically? I did not find something like:

QgsLayerTreeView.refesh()

in analogy to refreshing the map canvas. Am I missing something?


texts: Extracting Oracle geometry with FME's SQL Executor transformer

texts:
I am wondering why my SQL executor can't extract the geometry from my Oracle database. The geometry are polygons.

I added the 'SHAPE' column to the query in the executor, and I exposed the attribute but there is still 'no geometry available'.

Most of the time, this problem doesn't happen when using 'SQL Creator' (which in my case here is querying another dB).

How can I solve this?


texts: Line of Sight Analysis using Points and DEM

texts:
I have a DEM and a point feature class with the elevation value from the DEM as an attribute plus with x number of feet added to it representing a tower location.

I would like to do a line of sight analysis in all directions. Essentially, I want to create a buffer of x number of miles around the point layer, then using the DEM, perform a line of sight analysis using that buffer and the DEM to find areas in the DEM that are higher then the value in the point layer.

I have access to ArcMap/Pro with all extensions. I feel like the Observer Points tool is a good start, but it doesn't provide all of the functionality I need. Any ideas on how to accomplish this?


texts: Given distance (in meters) what is the angle in geographic reference systems?

texts:
In QGIS and other systems sometimes I have a layer in geographic reference system (typically WGS84, EPSG:4326) and I'd like to calculate a buffer of say 5000m around a point.

So, given a point coordinate (lat/lon) and a distance (in meters), what is the geographic length (angle) in a geographic reference system (for geographic lengths shorter than 180 degrees)?

I'm referring to two implementations (QGIS, postgis.net/docs/ST_Buffer.html ) which offer just that one "buffer radius" parameter.

Remarks:
 - I'm aware of Cosine Formula and Vincenty (like question 109562) but
   can't resolve the formulas to angle.
 - For the solution I prefer standalone implementations, say in Excel, pl/pgsql, Python or R.


texts:
I have Ubuntu 14.10.
I've installed osmosis from apt. I downloaded mapfile-writer from here

http://repo.typesafe.com/typesafe/simple/sonatype-snapshots-cache/org/mapsforge/mapsforge-map-writer/0.4.0-SNAPSHOT/mapsforge-map-writer-0.4.0-20140127.161606-1.jar

I put mapwriter to /usr/bin and /home/.openstreetmap/osmosis/plugins (just in case)

But when I run following command I get an error

/usr/bin/osmosis --read-xml file=UA.osm  --mapfile-writer file=dnipropetrovsk.map bbox=47.45,32.94,49.19,36.94 type=hd

:

SEVERE: Execution aborted. org.openstreetmap.osmosis.core.OsmosisRuntimeException: Task type mapfile-writer doesn't exist.
    at org.openstreetmap.osmosis.core.pipeline.common.TaskManagerFactoryRegister.getInstance(TaskManagerFactoryRegister.java:60)
    at org.openstreetmap.osmosis.core.pipeline.common.Pipeline.buildTasks(Pipeline.java:50)
    at org.openstreetmap.osmosis.core.pipeline.common.Pipeline.prepare(Pipeline.java:112)
    at org.openstreetmap.osmosis.core.Osmosis.run(Osmosis.java:86)
    at org.openstreetmap.osmosis.core.Osmosis.main(Osmosis.java:37)
    at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
    at sun.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:62)
    at sun.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
    at java.lang.reflect.Method.invoke(Method.java:483)
    at org.codehaus.plexus.classworlds.launcher.Launcher.launchStandard(Launcher.java:328)
    at org.codehaus.plexus.classworlds.launcher.Launcher.launch(Launcher.java:238)
    at org.codehaus.plexus.classworlds.launcher.Launcher.mainWithExitCode(Launcher.java:408)
    at org.codehaus.plexus.classworlds.launcher.Launcher.main(Launcher.java:351)
    at org.codehaus.classworlds.Launcher.main(Launcher.java:31)


texts: Excluding extent when polygonizing raster file using Python?

texts:
I want to polygonize a raster file. The raster file looks like this.

Here is the  raster file

Here is my codes in Python to polygonize the raster file. The codes source from Polygonize a Raster Band

#raster to polygons

from osgeo import gdal, ogr

#  get raster datasource

src_ds = gdal.Open( "C:\\Users\\select3.tif" )

srcband = src_ds.GetRasterBand(1)

#  create output datasource

dst_layername = "output_select3"
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource( dst_layername + ".shp" )
dst_layer = dst_ds.CreateLayer(dst_layername, srs = None )

gdal.Polygonize( srcband, None, dst_layer, -1, [], callback=None )
#The following two steps are very important, otherwise will generate error like (Coordinate not defined)
dst_ds.Destroy()
src_ds=None

I got the result which looks like this.

And here is the output shapefile

I do not satisfy with this, because I only want to have the polygon areas (Check mark) and exclude the extent (Cross mark) of the polygons. Could anybody tell me how to fix this in python?


texts: Getting SpatialPolygonDataFrame from gIntersection{rgeos] with sub-area as attribute in R?

texts:
I am using gIntersection (package rgeos) to intersect a polygon shapefile containing 20 buffers zones, with another polygon shapefile containing "ground occupation mode (GOM)" (an area divided in polygons where one polygon can be "forest", another "urban area", etc); both shapefiles loaded with readOGR, same projection.

What I want to get:

A SpatialPolygonDataFrame or Shapefile in which there is an attribute "ID "(corresponding to every buffer ID, therefore going from 1 to 20) and an attribute for every type of GOM in which the value is equal to the area of the said GOM for the corresponding buffer ID.

For a buffer area set to 450, the dataframe would look like

OR

What I actually get when doing gIntersection:

227 SpatialPolygon objects corresponding to each part of GOM polygon caught in the intersection process. The R console gives me a lot of "slots" informations such as area, coordinates etc. for each of these 227 polygons, so I guess the information is there, but I can't manage to make it a SpatialPolygonDataFrame or Shapefile with the informations displayed as previously shown.

I have no formation in spatial data analysis on R but what I learned by myself using Google, and I'm afraid this move is quite above my level, would any of you have a trick for me?

I tried your method and it didn't work quite well, let me show you the results:

plot(r)

text(r)

and plot(mos_cut,add=T) gave me the same black square  (mos_cut is the land use shapefile we used to call "gom")

I believe that my dataset is a bit different than the one you made up for the example, and this is my fault because in my previous post the dataframes and categories I displayed were a simplification of what my dataset really looks like. So I will show you what it really looks like:

Land use shapefile attribute table opened in QGIS, mos_2012_11 is the field for the land use category, there are 11 different categories (1= forest, 2= ...), and the dataset has around 20300 features for 3 fields.

Now this is the summary in R of the same shapefile when read:

This is the summary of my "buffers" layer:

If it can help, here's how I got the buffer layer:
(method applied to generate 20 random fixed distance buffers in my study area)

#reading shapefile containing the study area, basically a single buffer polygon
buffer_etude<-readOGR("C:/Users/JOLLYJUMPER/Desktop/R/Echantillonnage_R",layer="BUFFER1000saclay")
#transforming the shapefile into a point matrix to use the CSR function (generate random points in the polygon area)
buffer_etude<-slot(slot(slot(buffer_etude,"polygons")[[1]],"Polygons")[[1]],"coords")

colnames(buffer_etude)<-c("x","y")

buffer_etude<-as.points(buffer_etude)
#generating random points
random<-csr(buffer_etude,20)

colnames(random)<-c("x","y")
#transforming the random points layer into a shapefile to use the gBuffer function (to generate the said "buffers" layer we are talking about)
random<-as.data.frame(random)

random<-data.frame(Id=c(1:20),X=c(random$x),Y=c(random$y))

coords<-random[,c(3,2)]

attrib<-as.data.frame(random[,1])

names(attrib)<-"ID"

randomshapefile<-SpatialPointsDataFrame(coords, data = attrib, proj4string=CRS("+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs"))

#generating the "buffers" layer
buffers<-gBuffer(randomshapefile,byid=T,width=200,quadsegs=20,capStyle="ROUND")

then I basically apply your method (replacing "gom" with "mos_cut" and "buf" with "buffers" ) and I get these strange plots.
I tried to do the intersection anyway and after calculating for one or two minutes, R crashed to desktop.

I believe that my dataset structure is the root of the difficulties I'm having using your method...

Do you have any idea of what's wrong and how to solve this problem?


texts: Is there a way to customize the zoom button and legend in cartodb js?

texts: I am following this tutorial, http://docs.cartodb.com/tutorials/toggle_map_view.html, and have been successful in developing an application with similar functionality, currently I am working on customizing the appearance of the application. I would to customize the the appearance of the zoom button and the legend (font, color, location) in the application. Is there any wall to do this? In the css or js?


texts:
I've been trying to create visibility areas for CCTV in my city. Obviously simple geoprocessing tools aren't the answer, cause they create results as shown on the photo. CCTV can't look through buildings, so I was looking for some solution that will skip areas beyond building. Is there any tool or extension in ArcGIS Desktop?

It would be also very helpful to be able to create oriented buffers (since not all cameras work in 360degrees), but I haven't found any working solution yet.


texts: Analysing a network distance to two or more points using GRASS

texts:
I am doing a piece of analysis work looking at distances along a network in GRASS. Most of the analysis involves finding which parts of the network are within 600m of a series of points. GRASS does this very well using the v.net functions.

I've got stuck on trying to identify which parts of a network are within 600m of more than two points in the series.

Does anyone have any ideas how I could do this within GRASS?

I'm using the GRASS plugin in QGIS.


texts: Translating OSM 'osmc:symbol' into symbol geometry with QGIS

texts:
OpenStreetMap uses specific tags to describe route symbols that are used as waymarker (https://wiki.openstreetmap.org/wiki/Key:osmc:symbol).

For example: "osmc:symbol"=>"red:white:red_dot"

I would like to translate that key syntax into a symbol geometry with QGIS similar to the OSM composer does (https://wiki.openstreetmap.org/wiki/OSM_Composer).
But the only idea I get into my mind is:

build SVG-symbols for every possible combination (with Inkscape)
Name every symbol with the relating key (like 'red_white_red_dot.svg')
Use data defined override at an SVG marker line, where the SVG markers path is set up with the osmc:symbol field from the attribute table

That seems to be a not very elegant solution: Too static, too much work and too inflexible.

I'm looking for a better solution. Any ideas?


texts: Get a bounding box of a geometry that crosses the antimeridian using OGR

texts:
I'm trying to get a bounding box of the geometry that crosses the antimeridian using OGR:

from osgeo import ogr

geom = ogr.CreateGeometryFromJson('{"type":"MultiPolygon","coordinates":[[[[-180.0,-40.084015],[-180.0,-40.092945],[-179.988885,-40.092945],[-179.988885,-40.084015],[-180.0,-40.084015]]],[[[179.991288,-40.092945],[180.0,-40.092945],[180.0,-40.084015],[179.991288,-40.084015],[179.991288,-40.092945]]]]}')
minX, maxX, minY, maxY = geom.GetEnvelope()

The output is (-180.0, 180.0, -40.092945, -40.084015), which is not what I've expected. According to GeoJSON specification the bounding box should be (179.991288, -40.092945, -179.988885, -40.084015)

Is this a bug in OGR or my expectations are incorrect? Which is the best way to retrieve a bounding box in format defined in GeoJSON specification?

I'm writing code in Python using GDAL 2.3.1


texts: Earth Engine filtering - creating monthly composites, but returning UNIX timestamp rather than YMD

texts:
I'm new to Earth Engine and struggling with a problem related to ee.Date formats. I am filtering image collections to pass into later analyses, including regression. However, when I later try to add a time band, to permit that regression, I receive errors. We've tracked this down to the way we have defined system:time_start being wrong in the monthly filtering, because our code pulls the system:time_start in date format but really it should be the time stamp in milliseconds since the UNIX epoch.

Here is the code we are using to create our monthly composites:

var interval = 1;
var month_list = ee.List.sequence(1,12, interval);
var year_list = ee.List.sequence(2000, 2014);

var monthlyCHIRPS =  ee.ImageCollection.fromImages(
  year_list.map(function (ynz) {
  return month_list.map(function(mnz){
  var w = CHIRPSCollection.filter(ee.Filter.calendarRange(ynz, ynz, 'year'))
                        .filter(ee.Filter.calendarRange(mnz, ee.Number(mnz).add(interval), 'month'))
                        .sum();
  return w.set('year', ynz)
                        .set('month', mnz)
                        .set('date', ee.Date.fromYMD(ynz,mnz,1))
                        .set('system:time_start',ee.Date.fromYMD(ynz,mnz,1))
                        .set('ID',ee.Date.fromYMD(ynz,mnz,1))

})
}).flatten())

print(monthlyCHIRPS, 'Monthly CHIRPS');

I believe the error is in setting the system:time_start, in that currently we pull the date from YMD, in the same format. What I want to do, is somehow pull the date, but convert it back into a UNIX timestamp as in the original images, so I can later add a timeband.

Does anyone know how I can correct this code to do that, or have an alternative code for creating monthly composites that would permit this?


texts: Using a fixed map in the OpenLayers 3 Overview Map

texts:
Is there a way to 'lock' the overview map in OpenLayers 3 using a fixed extent?

The current implementation resizes the little rectangle and the map. So when I zoom in or out in my map, the map in the overview map might change as well.

However, I would like the map in the overview map to use a fixed extent, so that it never moves. So when I zoom or pan the map always stays the same and only the rectangle is adjusted.

The Overview Panel in QGIS provides a good example of what I mean.

In the following screenshots I am moving and zooming in QGIS, but my overview map always stays the same and solely the rectangle is adjusted.

I would like to implement that behavior in OpenLayers 3 or 4.


texts:
I'm using mapserver and openlayers for my gis server.
I have a sqlite database that contain cached tiles and i want to use this cached tiles as base layer.

databse structure :
    Tiles
        id
        x
        y
        z
        Type
        CacheTime
    TilesData
        id
        Tile(blob)

How can i do it?

Can i use my tiles as raster layer in MapServer? or as cache in MapCahce? or directly in OpenLayers?

I found MapCache, TileCache MapProxy and GeoWebCach(for openlayers). But i don't know Which can help me for this purpose.

note :
I can export my data to disk or change database tables, if needed.


texts:
I'm attempting to import gdal in Python 3.6 and I receive the error below.

  ImportError:
  dlopen(/anaconda3/lib/python3.6/site-packages/osgeo/_gdal.cpython-36m-darwin.so,
  2): Symbol not found: __ZN11xercesc_3_211InputSource11setEncodingEPKt
  Referenced from: /anaconda3/lib/libgdal.20.dylib   Expected in:
  /anaconda3/lib/libxerces-c-3.2.dylib  in
  /anaconda3/lib/libgdal.20.dylib

I currently use anaconda. I've attempted uninstalling/reinstalling gdal and libgdal with "conda install 'module'" with no success.

Conda is automatically installing version 2.3.0 for gdal and libgdal.


texts: How to graphically offset the boundary line of a polygon in QGIS?

texts:
In Arcgis, there is a way to offset graphically the contour line of a polygon like below. I'm trying to do the same thing in a styling method (i.e the solution can't be by buffering physically the layer). QGIS does have an offset option, but it's an XY type offset that moves the position of the line.

ANSWERS :

Way 1: use the geometry generator fill style. Just get in the box buffer($geometry,[your buffer distance]). The only difference is that Arcgis offsets the line whereas this solution creates a buffer, thus rounding the edges...

Way 2 : exact replica of Arcgis, inspired below by @iant answer : create in the style another style "Outline : Simple line". In this style, there is the offset box i was looking for ! It does exactly the same as in arcgis.


texts: Fill options for vector layers in leaflet

texts:
In leaflet, how can I manage the filling option? More specifically, I would like to fill inside a circle until a certain border not the whole geometry. Is that possible?
I simply use L.Circle: link

var annonciation_195 = L.circle([44.519, 6.760], 500,{color:'#C02900', weight:1, opacity:1, fillColor: '#D19D1F', fillOpacity:.5}).addTo(map);

I tried leaflet semicircle plugin link, but it also does not have this option to my understanding.


texts:
I want to assign projection to shapefile in my web C# application, for this i am using Gdal/OGR/OSR C# bindings and add osr_csharp.dll and ogr_csharp.dll references. But i am facing some difficulty while writing the code for assigning projection to shapefile. I want to know is there any API documentation for this which can help me to understand the classes/methods of these libraries.
OR if anybody have some piece of code for assigning projection to shapefile , then it will be great help for me.


texts:
How do I get the extent out of the GeoTiff format? I want something like this:

Extent(293518.1886150768,5680494.194041155,890338.5054657329,6267530.571271311)

AKA the xmin, xmax, ymin, ymax coordinates from the GeoTiff?

Edit: I'm implementing my own reader so I will need to know how the file format tells me this information.


texts: Downloading Census data to PostgreSQL PostGIS TIGER extension for Amazon RDS database

texts:
Can someone point me toward instructions on how to download census data for the postGIS_tigerP_geocoder extension?

I would like to use the geocode() function, but any query returns 0 results. From reading, I gathered this is because I must first download address and coordinate data from the US census website. Our database is Amazon RDS. My local machine's operating system is Windows 10.

I followed these instructions in the command prompt on my Windows machine ...
https://sparkgeo.com/blog/building-a-us-census-tracts-postgis-database/
When I got up to mget *, the ftp connection responded with

200 Switching to ASCII mode.
Cannot find list of remote files.


texts: Error mapping categorical vector from spatial points data frame in R

texts:
I have had trouble mapping categorical points on to a map using rworldmaps or the maps package in R.

The original data was Aragonite saturation levels in the ocean on an .adf file which was in Lambert Cylindrical Equal Area format ("+proj=cea +datum=wgs84 +lon_0=-160.0 +lat_ts=0.0 +x_0=0.0 +y_0=0.0"). When I plotted this data I get the below map.

However I wanted to convert these coordinates to lat/long and reclassify the Aragonite values into categories of 1 - 5. So I used the rasterToPoints(r, spatial = TRUE) function to create a SpatialPointsDataFrame that I reprojected to "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0" using the spTransform function.

head(SPDF2)
  #Aragonite     long      lat
  #1 1.542071 20.89833 84.66917
  #2 1.538187 22.69496 84.66917
  #3 1.537830 24.49159 84.66917
  #4 1.534834 26.28822 84.66917
  #5 1.534595 28.08485 84.66917
  #6 1.532505 29.88148 84.66917

At this point I reclassified the data onto my 1-5 scale by writing a CSV file, reclassifying the Aragonite values and reading back into R.

head(r.csv)
      long      lat Aragonite
1 20.89833 84.66917        5
2 22.69496 84.66917        5
3 24.49159 84.66917        5
4 26.28822 84.66917        5
5 28.08485 84.66917        5
6 29.88148 84.66917        5

I have created a SpatialPointsDataFrame

coordinates(r.csv) <- ~ long + lat

# And assigned the correct CRS
proj4string(r.csv) = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

#I then defined the vector I wanted to plot
plotr <- r.csv$Aragonite

# Defined number of colours to be used in plot
nclr <- 5 #the number of categories for Aragonite
mypalette <- brewer.pal(nclr, "Reds")

# Plot the map with desired lat/long coordinates and data points with colour coding
map("world", fill=TRUE, col="white", bg="lightblue", ylim=c(-60, 90), mar=c(0,0,0,0))
points(r.csv$long, r.csv$lat, pch = 15, col= mypalette, cex = 0.5)

HOWEVER, the resulting map is a mess! The points seem to be in the right place but the colours don't seem to match the new Aragonite categories at all?

I tried another tack and tried to rasterize my points but had no luck as they are irregularly spaced.


texts:
How do you appropriately setup a plugin data source to allow for null values (in .NET - C# preferred)?

I'm currently setting up a field using:

field.Type_2 = esriFieldType.esriFieldTypeSingle;
field.IsNullable_2 = true;
field.Name_2 = "Foo";

Then, later, during the cursor's QueryValues call, I'm trying to set a null value via:

// Given: public int QueryValues(IRowBuffer row)
row.Value[i] = null;

However, the value always maps through as 0, not null.  How do you properly setup the data source in a plugin data source to allow for a null value?


texts: How can I render the postgis layer to the openlayers map viewer?

texts: I have a postgres table which is storing the dynamic GPS point data. The geom column has been created. Is it possible to render the last 10 mins GPS data to the OpenLayers map? Is there an example?


texts: I am looking for a source of buildings elevation, preferably stored as 2D images (for each pixel, the intensity corresponds to the altitude in that XY). Considering the case of NY, it would be a XY projection (in greyscale) of this map. Do you know where I can find such material, in case it's available? And/or, is there any easy way to convert GIS data formats into images?


texts: SAGA Maximum Flow Length problems with the Coordinate systems

texts:
I'm new to SAGA.

I've got a srtm data set from europe and want to charge the maximum flow length especially from the river Elbe. The data is in WGS84
So first I ran the tool "Fill Sinks"(Wang & Liu) and the flow length downstream but the output flow length were in degree. So I transformed my original data in another coordinate system (I've tried EPSG 4087, 32662 and 54031) with meters as unit and ran the same tools but the watersheds and the flow length looks really creepy in those cases.

Also I have no clue what the "maximum slope legth" in the "Fill sinks" tool stands for.


texts: Raster reclassify using python, gdal and numpy

texts:
I would like to reclassify a raster file from a raster with 10 classes to a raster with 8 classes using Python, GDAL and/or numpy. The classes are represented as integers. I have tried following the steps from this post Reclassifying rasters using GDAL and Python? , the numpy.equal doc and also gdal_calc doc. However, to no avail.
The raster file to be reclassified has integer values ranging from 0 to 11 and also include values 100 and 255. The following show the reclass (from value : to value):
nodata : 4,
0 : 4,
1 : 1,
2 : 2,
3 : 3,
4 : 3,
5 : 4,
6 : 5,
7 : 5,
8 : 6,
9 : 7,
10 : 8,
100 : nodata,
255 : nodata,
What I have been able to do is select the raster file to be reclassified using tkinter.FileDialog and get the raster info such as geotransform, and pixel size with  reclass = gdal.Open(raster, GA_ReadOnly).
How do I go about solving the above.
It might be worth mentioning that the rasters to be reclassified can be fairly large in some cases (500mb to 5gb).


texts:
Michael Speth and i are in the process of fixing the PROJ.4 implementation of the HEALPix map projection (see http://trac.osgeo.org/proj/browser/trunk/proj/src/PJ_healpix.c for the current buggy implementation) and are running into some difficulties in the ellipsoidal case.

Here's the issue.  We've implemented the spherical version of the HEALPix forward projection via

FORWARD(s_healpix_forward); /* sphere */
    return healpix_sphere(lp, P);
}

and that works fine.  We now want to reuse healpix_sphere() in the ellipsoidal version of the HEALPix forward projection by simply replacing latitude lp.phi with authalic latitude and replacing the ellipsoid P with its authalic sphere.  What's the best way to do this in PROJ.4?

We tried

FORWARD(e_healpix_forward); /* ellipsoid */
    lp.phi = auth_lat(lp.phi, P->e, 0);
    P->a = auth_radius(P->a, P->e);
    P->e = 0;
    return healpix_sphere(lp, P);
}

where auth_lat() and auth_radius() are our homemade authalic latitude and radius functions. But the numbers coming out aren't correct.  We think that our attempt at changing P into its authalic sphere via the P->a and P->e assignments is failing.  Any tips?

Michael also asked this question on the PROJ.4 forum (http://lists.maptools.org/pipermail/proj/2012-November/006465.html) last week but hasn't received a response.

Thanks for your attention.


texts: Using OrfeoToolBox in Python 2.7.9 crashes python

texts:
I wanted to create a python script using Orfeo Toolbox for Orthorectification of AWiFS Data. I set the environment variables necessary and the following command runs fine from the command prompt:

otbApplicationLauncherCommandLine OrthoRectification -io.in "BAND2_RPC.tif" -map epsg -map.epsg.code 4326 -outputs.mode autosize -outputs.default 0 -elev.default 0 -interpolator bco -interpolator.bco.radius 2 -opt.ram 128 -opt.gridspacing 4 -io.out "BAND2_Ortho.tif"

But the following python script crashes when I try to run it:

import otbApplication
OrthoRectification =
OtbApplication.Registry.CreateApplication("OrthoRectification")
OrthoRectification.SetParameterString("io.in", "BAND2_RPC.tif")
OrthoRectification.SetParameterString("map","epsg")
OrthoRectification.SetParameterInt("map.epsg.code",4326)
OrthoRectification.SetParameterString("outputs.mode","autosize")
OrthoRectification.SetParameterInt("outputs.default",0)
OrthoRectification.SetParameterInt("elev.default", 0)
OrthoRectification.SetParameterString("interpolator","bco")
OrthoRectification.SetParameterInt("interpolator.bco.radius", 2)
OrthoRectification.SetParameterInt("opt.ram",128)
OrthoRectification.SetParameterInt("opt.gridspacing", 4)
OrthoRectification.SetParameterString("io.out", "BAND2_OrthoRect.tif")
OrthoRectification.ExecuteAndWriteOutput()

I am using OTB-6.4.0-Win64 and Python 2.7.9. The otbApplication module imports fine without any errors. When the ExecuteAndWriteOutput() method is to be executed, the program crashes with the dialogbox "python.exe stopped working".


texts: Raster calculator : How to make 0 + NoData = 0?

texts:
I have a point data layer (vector) that I converted to a raster layer. I want to replace the NoData value (between each raster point) in the point raster layer by 0. So, I made a background raster layer (value = 0).

My plan was to add the background layer of 0 to the point raster layer, but as we can see in the + (Addition) tool page, a number (0) + NoData = NoData.

Any tips on how to achieve this?

I'm using ArcGIS 10.3, but I also have access to QGIS 2.10.

EDIT : I saw the QGIS Raster Calculator: Need values to be set to 0 instead of no data question, but I can't find the "Reclassify grid cells" tool. I just have "Reclassify grid values".


texts: How does Empirical Bayesian Kriging of ArcGIS Pro handle sampling bias?

texts:
In ArcGIS Pro, simple kriging approaches allow you to de-cluster your data before applying a normal-score transformation. According to the Esri guideline, normal-score transformations are only feasible if the input data does not exhibit a strong sampling bias.
Empirical Bayesian Kriging does allow you to apply a normal-score transformation to the data before interpolation, but does not allow you to (at least manually) apply de-clustering algorithms on the data. Following that logic, EBK must either apply some sort of de-clustering technique automatically or must otherwise be suited to properly produce a prediction surface from geographically biased input data. Neither the Esri guidelines/publications nor other sources provide any information on this.
How then does EBK handle sampling bias?


texts:
I use the codes below to make an ajax call requesting points and show it on the map:

var layer_kos = new L.GeoJSON.AJAX("<?php echo base_url('assets/GIS/layer/request_kos.php');?>",{
pointToLayer: function (feature,latlng) {
    addRowTable(feature,latlng);
    tablesaw(feature);

    buffer(feature,latlng);//this function is used to add points to my own array
    return L.marker (latlng, {
        icon: L.AwesomeMarkers.icon({
            icon:'bed',
            markerColor: 'cadetblue' ,
            prefix:'fa'
        }),
        title: feature.properties.nama,
        riseOnHover: true
    });
},
filter:filterkos,
onEachFeature:function(feature, layer){
        layer.bindPopup(
            "
" + feature.properties.nama + "
"+
            "

 HARGA
 	  Rp. " +  feature.properties.harga + ",00" +
            "

 FASILITAS
 	  " + feature.properties.fasilitas + "

"
            // "<left> HARGA : Rp" + feature.properties.harga + ",00<left>
"+
            // "<left> FASILITAS : " + feature.properties.fasilitas + "<left>"
            )
}
}).addTo(map);

I need to works with the array obtained by buffer(feature, latlng), only when the ajax completed. I've tried on dataLoadComplete bellow but it doesn't seem to work.

layer_kos.on('dataLoadComplete',function(e){
        alert(buffers.length);
});

What is the callback function that I can use?


texts: Sync changes from view to feature class (using a GP tool; without fully replacing the data)

texts:
I have an ST_GEOMETRY view with 50,000 rows (in an Oracle 18c/10.7.1 geodatabase). The view references a WORKORDER table in an external system (via a dblink/Oracle 12c).
Problem:
The view is too slow to be used in ArcGIS; there is a 5-10 second refresh-time due to the dblink.

As an alternative to using the view directly (in real-time), I want to synchronize the view's resultset to a feature class (near-real-time).
I would sync any changes in the view to the feature class on a 1-minute schedule (scheduled task on the server).

I would only need to sync new or modified rows. Rows never get deleted in the external system — so I won't need to worry about deletes.
The changes to the data in the view are minimal: there is an average of only 5 rows changed per minute.
The logic would look like this: Sync rows from the view to the feature class where the changed_date in the view is newer than the last_edited_date in the feature class (or where the record doesn't exist in the feature class).

Question:
Is there a geoprocessing tool that I could use, on a schedule, to sync the data from the view to the feature class?

I think I would prefer the sync to happen at a row-by-row level, rather than doing a full truncate and replace (truncate & replace would seem excessive, especially on a 1-minute schedule).


texts: using openlayers 3 and three.js to use openlayers 3 map as texture

texts:
I am trying to use Openlayers 3 along with three.js to use openlayers output map as texture. I found that openlayers 3 is webgl based and so is three.js. So they must be compatible.

What I have done is displayed a map using ol.js function ol.Map().Now i want the map to be used as a texture in a terrain.

I think if i get the map as an image file instead of displaying in div container. The problem would be solved. But am unable to do it.

I have a terrain and a map but unable to use the map as texture for the terrain.

Can anybody help me with this. I am stuck and with no idea of what to do?
I have asked this in stackoverflow also.


texts: Time required to get results by solar radiation tool?

texts: How much does time does the ArcGIS solar radiation tool take to give solar map for day of about 50 km area?


texts: Moving and Rotating Raster Files in QGIS Individually?

texts:
I would like to overlap 2 raster files(.JP2/ECW) where images are captured from the same place but just with a little time difference. When I try to overlap the images by moving and rotating, QGIS rotates the whole canvas, so I could not move or rotate single raster layer. I tried the folllowing python code for this job:

from osgeo import gdal, osr

layer = iface.activeLayer()
provider = layer.dataProvider()
path = provider.dataSourceUri()
dataset = gdal.Open(path)

geotransform = dataset.GetGeoTransform()
gt = list(geotransform)
gt[0] = 500
gt[3]  = 500
gt[2] = 3 #x pixel rotation
gt[4] = 3 #y pixel rotation

dataset.SetGeoTransform(tuple(gt))

dataset = None

but it does not update the current map on the canvas.


texts: Merging 4 non-referenced images using QGIS

texts:
I have four .tif images of the same (MxN) size, which I want to merge together. The images are not georeferenced (i.e., pixels are not associated to latitude/longitude coordinates). The 4 images are part of a same global map of the Mars surface and correspond to the top-left, top-right, bottom-left and bottom-right "portions" of such global map.
When I load the 4 images as raster layers in QGIS, they all overlap, since they are not georeferences. Therefore, when I use the function Merge, I obtain a (MxN) image instead of the (2Mx2N) desired image.
Is there a way to easily merge them in the correct way (i.e., telling that image 1 is the top-left one, image 2 is the top-right one, etc.)?


texts: How to import a KMZ file and load data into correct columns in QGIS

texts: I have downloaded building heights attribute data from Ordnance Survey Mastermap. The format of the file is KMZ. After importing into QGIS, the data I need was displayed in one column named as "description" rather than individual columns as shown in the figure below. How can I split the description column into several columns to get the data I need?


texts: How to use dissolve a subset of a PostGIS table based on a value in a column

texts: I have a spatial table in PostGIS that has millions of land parcels. The parcels have ownership information and also a column that specifies whether the parcel is private or public. I would like to merge/dissolve only the private parcels that have the same owner name and share a boundary (adjacent parcels owned by same owner) but also keep the private parcels that were not dissolved.Below is a sample table. Based on this table I would like to dissolve parcels with id 2,3, and 4 into a single parcel (assuming they are adjacent) but also keep parcel with id 1. I tired using ST_Union but when I specify the union on column owner_name, all the parcels (including parcel id =1) just have one column owner_name in the final table. I would like to keep all the columns for the parcels that did not merge/dissolve and keep only owner_name for the merged/dissolved parcel. So parcel with id 1 should have all the columns and the merged parcel from 2,3,4 should have only owner_name. How can I achieve this? Ideally I will like to keep all the resultant data in one table.


texts: Joining attributes between layers based on proximity of features

texts:
In QGIS 3.12 I have a Layer Y that has a specific attribute and my objective is to fill this attribute with information from Layer X regarding the proximity between the features of Layer X and Layer Y.
I'll exemplify below:

Point A and B are from Layer Y and Point 1,2,3 and 4 are from Layer X. I would like to know if it is possible to fill the attribute 'Address' from Layer Y with the attribute 'Address' from Layer X regarding the proximity of the features.
In this particular case I want to associate the 'Address' from Point 4 to the 'Address' from Point A, and the 'Address' from Point 3 to the 'Address' from Point B using the Field Calculator.


texts: Customizing OpenStreetMap "carto" style?

texts:
I have been following the switch2osm tutorial explaining how to setup your own tile server:

https://switch2osm.org/manually-building-a-tile-server-16-04-2-lts

Everything works fine, however I would like to tweak a few things in order to better understand the rendering curtomization process. As a start I would like to change the color of the country bounderies from the default openstreetmap color to black. Also the bounderies are only displayed at low zoom levels (such as 3) but I would like to enable them at any zoom level.

My understanding is that this kind of styling is controlled by the Mapnik.xml file, which is in turn generated by the carto project.mml > mapnik.xml command.

How do I modify the project.mml or one of its dependent files in order to achieve my goals?


texts: Convert maps from different coordinates scale in Cartesian coordinate system

texts:
I need help to convert position of objects on map to a new scale.

Coordinate (0, 0) is the top left corner of the map.
Currently, 1 point distance is X
I need to convert 1 point distance to be Y
I need to update the location of each object on the map to the new scale

What I was thinking to do is this:

The scale is X:Y, which is basically x/y. Let K = X/Y
Fore each object on the map with coordinates obj_x and obj_y, new coordinates are calculated the following way:

new_obj_x := obj_x / K

new_obj_y := obj_y / K

Is this correct? I have a feeling that something is very worng


texts: How to use onEachFeature in Leaflet.js map?

texts:
I have a JSFiddle and I am trying to implement onEachFeature (at the bottom) of the code.

Each polygon has its own popupContent as in the onEachFeature tutorial shown on Leaflet.js.com.

However, I get stuck trying to merge together these pieces of code:

var myGeoJson = L.geoJson(myData, {
  style: style
}).addTo(map).bindPopup("test");
var popup = L.popup();

function onEachFeature(feature, layer) {
    if (feature.properties && feature.properties.popupContent) {
        layer.bindPopup(feature.properties.popupContent);
    }
}

var myGeoJson = L.geoJson(geojsonFeature, {
    onEachFeature: onEachFeature
}).addTo(map);

Ultimately, I'm just trying to pass in the features.properties.popupContent as a variable any way possible without losing my style properties (fill color, etc).


texts: How to tell if an ArcGIS Service is secure?

texts:
I am creating a .NET REST SOE for ArcGIS Server 10, and I will need to determine if a map service (or image) is secure or not.

Is there a property on a MapService that would tell me this in either ArcObjects or REST?


texts: Set a RANGE of values as nodata in a .vrt?

texts:
Is it possible to set a range of values to nodata in a vrt file? How?

I have some elevation models with incorrect values, down in the -15,000s. I could just clip the rasters along that edge, or apply some raster calculation to transform them to nodata, but that's a fair bit of monkeying around. It'd be nicer to just say "everything less than -100 is nodata" and not duplicate the 99.9% of the cells that are just fine to fix the bad ones.

In reading the GDAL Virtual Format Tutorial I see a ComplexSource option to add "a custom lookup table to transform the source values to the destination. ... following form: <LUT>[src value 1]:[dest value 1],[src value 2]:[dest value 2],...</LUT>

but when I try that it transforms all values, not just the ones within the range.

Min/max from source: Computed Min/Max=-16111.000,2618.000
VRT Lookup transform: <LUT>-32767:-10,-1000:-1</LUT>
Out min/max: Computed Min/Max=-5.000,-1.000

Am I doing it wrong, or is what I'm attempting simply not supported? (by any VRT construct, not necessarily just LUT). Is there some other method of setting a range to nodata that doesn't require duplicating files to get there?

Source data:http://www.viewfinderpanoramas.org/dem3/P07.zip (34MB, home page), any tile ending in W142.