texts: Match SpatialPoints to Delaunay Triangles in R

texts:
I have 2 SpatialPointsDataFrame files. I compute the Delaunay Triangles from the first shapefile (shp1) as in the code below. Then I wish to assign each point from my second shapefile (shp2) to one of the Delaunay Triangles held within the object vtess.

vtess <- deldir(shp1$coords.x1,shp1$coords.x2)

I can assign each shp2 point to a tessellation tile in the following way:

tl <- tile.list(vtess)
shp2$tile <- 99
for (i in 1:length(shp2$tile)){
   shp2$tile[i] <- which.tile(shp2$coords.x1[i], shp2$coords.x2[i], tl)
}

However, while I can find the list of triangles using triang.list(vtess), there is no which.triang command, so it is not possible to similarly assign shp2 points to triangles.

I'm guessing I therefore need to (1) create a SpatialPolygon of the triangles, and then (2) overlap the shp2 points to the new polygon file. Is this the most efficient way? And if so, how do I create a polygon shapefile from the triangles held within vtess?

Cross-posted, with no reply, here: https://stackoverflow.com/questions/51700872/match-spatialpoints-to-delaunay-triangles


texts: Using ModelBuilder for adding and calculating field

texts:
How can I add a field, calculate the area in it then export the maps as an excel sheet using the ModelBuilder for all the data stored in the geodatabase file?
Note, the picture was my try to do the task.


texts: Add basemaps in QGIS - what differs Tile Services from QuickMapServices QMS?

texts:
Environment
QGIS-version: 3.14.0-Pi
QuickMapServices Plugin: 0.19.11.1

I want to use background maps inside QGIS (QGIS version 3.14.0-Pi) and can do that by using f.e.:

WebMapTileService (can be done inside QGIS as defined XYZ Tile Services)
QuickMapServices (are integrated inside the QMS Plugin)

Often as I tried to use QMS integration the first thing I see is a watermark or the hint to an expired API key. Sure I can use the catalog search and find many entries, but it looks like the catalog data are coming from nextgis.com (and they are expecting a key which isn't possible for me to enter somewhere inside QMS).
Whereas as I'm using the WebMapTileService nothing similar has happened - and I'm also able to pass a key into the request.
So I'm a bit confused and don't understand:
What is technically the difference in between QMS versus WebapTileService?


texts: How to convert Layout View to DXF or DWG?

texts:
How can I convert Layout View to a DXF or DWG format?

I have tried the PDF way i.e.

Export to PDF
Convert from PDF to DWG using PDF2cad

It was not a great experience because the DWG does not seem good and is not georeferenced.

Can anyone offer help or experiences to help solve this issue?


texts: Change zoom factor of QGIS programmatically with PyQGIS

texts:
I would like to change the zoom factor of QGIS dynamically without restarting QGIS.
The changed settings I would like to use for a plugin would be that it let's me dynamically change the zoom behaviour when I use the mouse wheel for zooming in and out.

I know that I can change the QGIS-Settings like that:

from PyQt4.QtCore import QSettings
s = QSettings() #getting qgis settings
s.setValue("/Qgis/zoom_factor", 1.1)
s.sync()

But this requires a restart of QGIS to really change the zoom behaviour.
Is there anything else I could do to make QGIS accept the new settings like it does when I click the "OK"-button in the settings dialog?


texts: Calculating Mean for Multiple Rows in ArcGIS for Desktop?

texts:
I am trying to calculate a mean across a multiple number of attribute rows and have this value added to a new field that I have entered in the attribute table. I have an attribute table which has a number of different ecoregion types associated with them.

For example:

If I select "Arid Canyonlands" there are 3 attribute rows that are selected, each with their own "mean". I am trying to add a new field to the attribute table which will take the "mean" of these three and populate that value in the new field created. I could manually do it through an editing session and using a calculator but just wanted to know if there is an easier way to actually do it through field calculator?


texts: Implementing layer.setEditType for custom widgets in PyQGIS?

texts:
How could I use layer.setEditType() method in QGIS for a custom form widget? Cusom form elements are well explained by Matthias Kuhn here.

I was able to change the existing default widget types. For instance,
layer.setEditType(idx, layer.ValueMap) but when I tried the same on a custom widget, it doesn't work. How could I do it?


texts: Iterate through list of field names and calculate fields with PyQGIS?

texts:
I am using QGIS 3.2.1 Bonn 64 bit on Windows 10. I would like to loop through a list of field names and perform calculations on each field. Only the first loop in the iteration performs any calculations. The rest are skipped.

result_path = r"D:\QGIS\hex_results.geojson"

tracks = QgsProject.instance().mapLayersByName("patrol_tracks")[0]
result = QgsProject.instance().mapLayersByName("hex_results")[0]

t = tracks.getFeatures()
r = result.getFeatures()

r_prov = result.dataProvider()
r_fields = r_prov.fieldNameMap()

months = [x for x in processing.run("qgis:listuniquevalues", {'INPUT':tracks,'FIELDS':"Month_Year",'OUTPUT':'u'})['UNIQUE_VALUES'].split(";")]

for m in months:
    result.startEditing()
    field_index = result.fields().lookupField(m)
    print (m)
    print (field_index)
    a = {field_index : 0}
    for feat in r:
        fid = feat.id()
        r_prov.changeAttributeValues({fid : a})

    result.commitChanges()

The print calls work for every list item but the field attributes are only altered for the first iteration.

Another method I have tried (using the same variable names as above) is:

for m in months:
    processing.run('qgis:fieldcalculator', {'INPUT':result, 'FIELD_NAME':m, 'FIELD_TYPE':1, 'FIELD_LENGTH':10, 'FIELD_PRECISION':0, 'NEW_FIELD': False, 'FORMULA':'0', 'OUTPUT':result_path})

But with the same result - only the first field in the list is calculated. This code will be part of a more complex script that calculates field values based on the intersection of selected features - which is more or less working up until this point.

Any suggestions? I wondered if the dataProvider needs to be refreshed or updated after each iteration but commitChanges(), added later, doesn't seem to have any effect.


texts: Give different shapes the same color based on attributes

texts: I am working with ArcMap 10.4 and have some issues. I'm making a Floodmap based for Europe and want to color the shapes related to the total amount of floods in every polygon. For example I have Czech Republic where 3 of the 8 shapes has 5 floods as a total amount in the attribute table. Is it possible to give these 3 shapes (and all the other ones with 5 floods) automatically the same (orange) color?


texts:
I  have had no luck getting the buffer tool to accept anything but degrees as units of measure.

I have found lots of stuff saying the layer needs to be reprojected and saved but it hasn't worked at all for me.

Is there a way I could create a buffer without using ftools or at least force the units to meters somehow?

As a workaround I converted meters to degrees (lat) and used that but the final product needs to be as close to reality as possible.

Things I've tried:

setting every unit option I could find to meters (where possible).
setting everything to NAD83/Maryland (data is for Washington, DC) and saving it as such (as layers in ESRI shape files).
reimporting the reprojected layers
setting relevant layers to Google Mercator

The was tried followed by creating a buffer. Many were tried in combination.
QGIS 1.7.3
Slackware64 current (qgis from SBo-13.37 repo, tried on multilib and plain 64it with same results)


texts: How to use nextPageToken to check all files under one folder with more than 1000 files in Earth Engine?

texts: Recently I carried land cover classification and export around 3000 files into one folder of the asset. However there were only 1000 files displayed on the code editor. I also used ee.data.listAssets() function and print the result. The result showed a list of 1000 elements with nextPageToken code. I wonder how to use the nextPageToken code to check the rest file ids but I did not find the resolution. Other methods which could work are also welcome.


texts: Web Map is not displayed in Openlayers with no errors

texts:
I want to display a map from a geojson file. Below are my codes. The map didn't show and there were no errors in the console.

  var urlModev ='maps/StaRita.geojson';

  var styles = [
    new ol.style.Style({
      stroke: new ol.style.Stroke({
        color: 'blue',
        width: 3
      }),
      fill: new ol.style.Fill({
        color: 'rgba(0, 0, 255, 0.1)'
      })
    }),
    new ol.style.Style({
      image: new ol.style.Circle({
        radius: 5,
        fill: new ol.style.Fill({
          color: 'orange'
        })
      }),
      geometry: function(feature) {
        // return the coordinates of the first ring of the polygon
        var coordinates = feature.getGeometry().getCoordinates()[0];
        return new ol.geom.MultiPoint(coordinates);
      }
    })
  ];

        var vectorLayer = new ol.layer.Vector({
            source: new ol.source.Vector({
                url: urlModev,
                format: new ol.format.GeoJSON(),
                style: styles
            })
        });

        var EPSG ='EPSG:4326';

        var projection = ol.proj.get(EPSG);

              var map = new ol.Map({
                layers: [vectorLayer
                ],
                target: 'map',
            renderer: 'canvas',
                controls: ol.control.defaults({
                  attributionOptions: /** @type {olx.control.AttributionOptions} */ ({
                    collapsible: false
                  })
                }),
                view: new ol.View({
                projection:projection,
                  center: [19948.84, 19955.93],
                  zoom: 1
                })
              });


texts:
I have already built Ubuntu for OpenStreetMap and I have some LiDAR data file in LAS/LAZ format. Is it possible to make LiDAR data appear in 2D and be a part of OpenStreetMap? If it is, how?
I have already pushed LiDAR to PostgreSQL pointcloud. Make surface model by lopocs and Cesium
And I try some tools to handle LiDAR like liblas, lastool, lidar2tiles,...


texts: Connecting points located in opposite side of international date meridian using Python GDAL/OGR?

texts:
Given a list of coordinates, I need to create a line shapefile joining each of the point with a Python GDAL/OGR script.
The procedure described here Convert XY points to a line?
 works fine in almost all cases, except when the points are crossing the date change meridian. When this happens, the line connects the two points by crossing the entire world, which is a obvious mistake.

How can I force the script to connect the 2 points through the shortest distance?

Here are some coordinate examples:

lon, lat
172.5, 50.0
175.0, 43.4
-175.0, 35.0
-173.1, 31
-172.1, 28.0
-171.3, 26.5
-173.9, 21.8


texts: FME - How to send list of files to reader

texts:
I have an Excel spreadsheet which contains the file paths of MapInfo files that I want to convert to a File Geodatabase using FME.  Is it possible to send this spreadsheet to the MapInfo reader so that it only converts those files?

So the spreadsheet might look something like this:

File Path
C:\Temp\File1.TAB
C:\Temp\File2.TAB
D:\Folder\File3.TAB
G:\Other Folder\File4.TAB

All file paths are in column A of the spreadsheet.


texts: Populating Fields in Table with GUIDs from Feature Class?

texts:
Is there a way to add GUIDs from a feature class to a table based on specific attributes? I have a large maintenance record table that is tied to specific points through the GUID of that point. I have added new records but they are not tied to a point yet. Instead of going through and copy/pasting each pole GUID to the appropriate maintenance record manually, I was hoping to use model builder or python to populate the fields in the table automatically. The data that overlaps between the two is the name of the line and the pole number.

I think this is a relationship (many to many), but both the line name and the pole number need to match. I'm also new to using model builder and python.


texts:
ArcSDE 10 with a MS SQL 2005 backend.

I am trying to see the most recently modified geometries using SQL so that I can use it to update another GIS database running PostGIS. Has anyone tried this before? I am using jTDS and oracle SQL developer to connect to the MS SQL database, if that helps.


texts: Is it possible to simplify orthogonal polygons with OpenCarto Java Library?

texts: I want to know if simplification of orthogonal polygons can be done with OpenCarto Java Library. If yes, does it have any API documentation?


texts:
When I export a svg file I get some wrong rendering, like maps exceeding their frame, grid lines that are supposed to be disabled, and inconsistent line styles.
This happens only when exporting in vector, rasters are perfect as seen in the composer preview.

What could cause this?


texts: Formatting OSM Data For Use in ESRI Address Locator

texts:
I've been given the task of researching the possibility of setting up an ArcGIS Geocoding Service that can work with OSM data.

I am aware of Nominatim, but since these service(s) will be running on Windows Server 2008, and the customer has already invested big $$ on the ArcGIS Server license - I don't think its going to be an option.  From what I've read, installing Nominatim on Windows is a daunting task in itself.  If someone knows of any tutorials of how to accomplish this, please let me know.

I am far from an expert in either geocoding or OSM, but I would think the biggest barrier to using OSM data in an ESRI address locator is lack of address ranges.

Is there any way at all to get address ranges for a given OSM road segment?  Any other tips or suggestions welcome.


texts: Reprojecting GeoTIFF from WGS84 to UTM with gdalwarp?

texts:
I have a GeoTIFF image in WGS84 and when I try with gdalwarp to reproject to UTM, it has a problem.

Size of TIFF is  J: 15769, W: 60027

J_LeftTop：2.0705046510  W_LeftTop：0.4974727391；
J_RightTop：2.0702337994 W_RightTop：0.4996621899；
J_LeftBottom：2.0760866619    W_LeftBottom：0.4980068600；
J_RightBottom：2.0758219229 W_RightBottom：0.5002002912；
J_CNT：2.0731617464 W_CNT：0.4988380007

When I try:

gdalwarp -t_srs '+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs' test.tif test_utm1.tif

  ERROR 1: Too many points (10201 out of 10201) failed to transform,
  unable to compute output bounds.

What should I do differently?


texts: QGIS Print Composer: Make grid lines around single layer?

texts:
Using QGIS 2.12, I´d like to know if it´s possible to make grid lines only around a certain map layer instead of the whole worksheet?
I prefer having legends and scale bars outside the map, which makes it annoying when the grid lines and labels continue outside of the map layer.


texts: Shifting points along normal vector to path

texts:
I have list of points which are connected with line in "zig-zag" fashion:

Using QGIS, I'm trying to create a line connecting all the middle or inside points and then shift all outside points to this middle line along normal vector.
End result is a single line with all points plotted on this line.

Note that intervals could be irregular and I can't simply divide middle line in equal portions and create new points.

Which QGIS tools would be appropriate for this task?


texts: How to convert a shapefile from bessel1841 to WGS84?

texts:
I would like to convert a shapefile which has +ellps=bessel to +ellps=WGS. In my understanding, even QGIS 1.8.0, it is impossible.

Let me know the free software which is able to convert shapefile from bessel to WGS.


texts: delete unconnected polygons in a shapefile

texts: I need to delete isolated polygon (green circle) from a shapefile but I do not want to delete polygon that connected with other polygon (red circle?


texts: Manually entering attribute information for each point in QGIS edit session?

texts:
I am new to QGIS and am editing an existing point shapefile that has 5 TYPE classifications. In ARCMap I can click on a symbolized TYPE point and the new digitized point will have the same attributes. In QGIS it seems that each newly added point in the file requires a manually attribute information.

Is there a tool or set-up that mimics the ease of editing in ArcMap?


texts: Creating simple polyline geometry objects with ArcPy?

texts:
I'm trying to learn how to create a polyline object from a set of XY pairs. I need to create an empty point object and an empty array object, place the coordinates in a list, and use a for loop to iterate through the list and populate the point object. From there I'm pretty sure I can create the polyline from the array, which is the last step.

How do I go from the list of points to the point object using a for loop? I've been trying to use the .add or .append functions to no avail.

# Import module
import arcpy

# Create empty objects:
#
# Point
point = arcpy.Point()
# Array
array = arcpy.Array()

# Create a list of coordinates (X, Y)
coordinates = [[260, 320], [450, 420], [360, 680], [600, 1000], [900, 190]]

# For loop iteration through coordinates
for pair in coordinates:
  point.x = pair[0]
  point.y = pair[1]

  # What else goes here to populate the point object?

print point # to test if loop is functioning correctly


texts: I am trying to create an NDVI image using an RSGISLib python script, but it is displaying no data values

texts:
My script runs fine and creates an output, but the output is displaying values all the way down to -999, rather than just having the values present in the NDVI result (0-1). I can create an NDVI image from this, but it requires some tweaking within a GIS program, which i'd rather not have to do as i have to process a large number of images. Is there anything i can add to the script so this stops happening?

This is the script

import rsgislib
import rsgislib.imagecalc
import rsgislib.imagecalc.calcindices

InImg = './LS5TM_19890814_lat55lon10028_r22p34_vmsk_rad_srefdem.kea'
OutImg = './Manitoba1989'
ndviImg = OutImg+ '_ndvi.tif'

rsgislib.imagecalc.calcindices.calcNDVI(InImg, 3, 4, ndviImg, stats=True, gdalFormat='GTIFF')


texts: Getting Error 000212: Cannot create XY event source Failed to execute (Make XY Event Layer). using ArcGIS ModelBuilder

texts:
I have a model created in ModelBuilder using ArcGIS Desktop that works on PC1 set up with relative paths and workspaces. When I try run it on PC2 I get an error message statin

  Error 000212: Cannot create XY event source Failed to execute (Make XY
  Event Layer).

with the exact same files. I tried again an run it on P1 again with the same files and it works.

Can you advise?


texts: Buffer in QGIS/R is 650-670 m wide instead of 500 m

texts:
I am trying to make a 500 m buffer for my polyline layer. Unexpectedly, I receive buffer varying from 650 to 670 m, depending on location. As many resources recommend to use UTM before making a buffer, I use UTM projection for zone 31N - EPSG:23031 - ED50 / UTM zone 31N (actually, the layer is already provided in this CRS).
Whole polyline layer is located in Catalunya (SE Spain). Even the buffers placed near in the very SE of Spain, which is right in the middle of the zone 31N (where no distortion should happen), are 650 m!

PS: I did this in QGIS and also in my script in R, receiving exactly the same result in both cases.
EDIT: using QGIS measure line tool (View > Measure > Measure line)


texts: Automatically Show the Page Number of Atlas QGIS

texts:
How can i show the page number of a atlas on QGIS?

I don't want to use an attribute, i just want to QGIS to automatically show de page number.

I believe that is a simple tool, but i Can find a way to do that.


texts: How to overcome field truncation when converting OSM to shapefile with QGIS?

texts: I am trying to convert an OpenStreetMap file into a shapefile, but I am having problems with the "tags" field being cut-off after 80 characters. With QGIS 1.4, I use the OpenStreetMap plugin to 'Load OSM from file'. The file loads in QGIS. If I right-click on the 'lines' layer and select 'Save as shapefile', I end up with a shapefile but one that has an 80 char limit on "tags". How can I increase this 80 character limit?


texts: point thinning/point selection based on distance threshold in QGIS/GRASS

texts:
I need to thin a point vector later based on the distances amongst points. So I have a threshold and wish to omit all records but one that is below a distance D representing a distance between points. Say D points are tree locations D is the canopy radius and I want to select/extract/remove all trees within the canopy of another tree.

I have tried 2 approaches unsuccessfully:

Distance matrix (QGIS) and I can select points below D but can only remove all or none, but I need to have a single remaining one.
Buffering the points with D and applying the Select Within Plugin to select polygons with their centroid within other polygons returns a selection of all polygon because centroids fall inside the other polygons.

How do I remove all records but the first that are closer than D?


texts: Match Button Greyed Out in Interactive Rematch Window

texts: I used ArcGIS 10.3.1 and built a composite address locator from 5 address locators. All locators are in the same file gdb. I geocoded addresses in a table within the same file gdb and tried to rematch unmatched results. For some of the results there is a 100 score candidate after I fixed the address string. But the match button (as well as the rematch button) was greyed out so I am not able to rematch this address with the candidate. A screenshot of my rematch window is as below:


texts:
This has been probably asked many times before; looking for answers to no avail. After first click this is what I get ...

then this ...

I'm still stuck with QGIS 2.18 'Las Palmas'

I'm not savvy with programming or computer science or anything in that sphere. I'm a landscape architect who happens to use QGIS for variety of analyses but definitely not as intense as GIS professionals.


texts: Running Python 2.6 and 2.7 in ArcGIS Desktop?

texts:
I figure that there must be a way for me to run python 2.6 and 2.7 on the same computer.  I've got a program that requires 2.7, but ArcGIS 10.0 requires 2.6.  I'd like to be able to run both programs at the same time.

I've installed both and in the geoprocessing settings set the default editor to C:\Python26\ArcGIS10.0\Lib\idlelib.idlepyw.  Now when I right click in ArcToolbox, select edit, the IDLE 2.6 editor comes up but the script won't appear like it formerly did.

I can view the script in IDLE 2.7 but that doesn't do me any good.


texts: GDAL with PostgreSql LIST_ALL_TABLES fails?

texts:
I have compiled GDAL 2.0.0 and am accessing a PG database. I note that the online documentation has two options available with the word "TABLES" in it. I think the logic in the code isn't handling the second syntax correctly.

(1) TABLES=string: Restricted set of tables to list (comma separated).

(2) LIST_ALL_TABLES=YES/NO: This may be "YES" to force all tables,
   including non-spatial ones, to be listed.

When I use a connect string like this:

PG:dbname='OpenGIS' host='xxxx' port='5432' user='xxxx' password='xxxx' active_schema=public schemas=public,tiger LIST_ALL_TABLES=YES

I get an error (my own handler in CPP program):

[PQconnectdb failed. missing "=" after "LIST_ALL_" in connection info string]

What I see in ogrpgdatasource.cpp at line 440 is some code starting like this:

pszTableStart = strstr(pszConnectionName, "tables=");
if (pszTableStart == NULL)
    pszTableStart = strstr(pszConnectionName, "TABLES=");

This code strips out the TABLES=... options for some reason, but the logic seems to ignore the other possibility of LIST_ALL_TABLES=YES.

I tried using mixed case options but that gives another error:

[invalid connection option "LIST_all_TabLeS"]

I've tried single quotes, double quotes etc but the result is the same.

What I am suspecting is that the parser is pretty simple and is not quite handling this right.

There may be a similar issue in SpatiaLite driver ? - I am working on that next.

I am working on Windows 10 and VS 2013 and don't have a project to build and test with.

Does anybody have a similar issue or am just using the wrong syntax?


texts: How to import OS street view data into MapInfo Pro

texts:
I'm trying to import multiple files into MapInfo Pro to create one tab file which is a map background of my local authority.

I'm not a regular user of MapInfo and my skills are at a beginners level. The files downloaded from Ordnance Survey contain Tab files in a georeferencing folder, jpegs in a data folder and tiff files in a data folder.

I have moved the tab files and jpegs into one folder together and I am able to view my map background. The layers are all fitting together and are geographically where they should be.

At this point, the layers are not editable or selectable. I did try saving a copy and opening that tab file but it remains unselectable and uneditable which means many of the menu options are unavailable to me.

I would like to know how I can combine all of these layers to make one tab file?


texts: Combining two raster layers into one vector points layer with two attributes

texts:
I have two raster layers, which I need to use to make a scatterplot. I can't find any tools that will let me produce a scatterplot from two raster layers, but there appear to be a few tools for producing a scatterplot from vector data.

So I'm trying to convert the two rasters into one vector points layer, with the values from each raster stored as attributes against each point. Is this possible, and how would I go about it? I've been googling all afternoon and can't work it out.


texts: Floodplain Delineation using Rasters - ArcGIS/HecRAS

texts:
I exported my HecRAS results as GIS format and created from the .sdf file a .xml file with the HecRas extension for ArcGIS.

After that I want to perform

"Floodplain Delineation using Rasters"

with the HecRAS extension.

But I got always the same error and don't know how to handle that.


texts: Filling Field2 Based on Values in Field 1 using Python Parser of ArcGIS Field Calculator?

texts:
I have been trying to rework attributes in a polygon shapefile and think I have figured out how to ask in a sensible manor.

I need to fill column Field2 with a 1 IF the value in column Field1 is equal to a certain value.

For instance I want a 1 in my field "Greenspace" if my field "leisure" has the value 'park' or 'golf_course' in the column.


texts: Buffering one side of line using GeoPandas

texts:
I have these transects I've buffered as training data, but now we had "edges" defined where we need a simple buffer on ONE side of the line. I run my scripts within Jupyter Notebooks, so it needs to work standalone. Right now my buffer code looks like this:
lines['geometry'] = lines.buffer(1, cap_style=3)
lines.to_file(linesBufferFolder + "\\Lines_Buffer")
This creates a 1m buffer on both sides with end caps. Like this:

I essentially need to bisect that polygon, and only keep the "left" side based on the order which the points go. My initial thought was to use a one sided buffer, but I can't find any documentation for that in GeoPandas Geometric Manipulation documentations.
My other idea was to simply erase half of the polygon I currently have, but again, to get that overlay I would think I would need to buffer the other side, bringing me back to my original issue. I know how to manually do one side buffers, but finding a Python library that uses it other than ArcPy is stumping me.


texts: Algorithm behind custom tool of ArcGIS for Desktop called xyz2grid21?

texts:
I am a researcher from Romania and I have used  xyz2grid21 ArcGIS tool (by Min-Lang Huang) in order to convert xyz files (LIDAR data) to ASCII and ESRI grid in order to obtain a DTM.

I need to know what algorithm does the XYZ2GRID tool use in order to rasterize the data? Is it using a mean or minimum value or does it interpolate somehow?


texts:
Is there an available patch for importing GeoPDFs, like the topographic maps that are available through USGS, into a session?

Software that costs $3000 is not a viable option.


texts: Reprojecting from Clarke_1880_RGS_Lambert_Conformal_Conic to ETRS_1989_LAEA?

texts:
I am unable to reproject Marseille shapefiles from Clarke 1880 RGS_Lambert_Conformal_Conic to ETRS_1989_LAEA.

I know that Clarke is a geographic coordinate system and it has no datum.  I have tried varying different things (project to ED50 utm zone 31N, ED79, WGS84) but it seems that ArcGIS does not do anything ... the data always remains in the same place.

The funny thing is that I have one shapefile in GCS_Clarke_1880_IGN that projected well to ETRS89 on the fly.


texts: Exporting layout elements and constants to Python script?

texts:
Is it possible to export the elements and constants used in a layout to a Python script? Something similar to creating a model and exporting it to Python but instead of being a model, it's a layout?

I have created a layout which contains 4 dataframes.  I would like to create a python script which loops through folders and updates the layout based on the layers in each folder (each folder represent a species). I also want each iteration to export to a .png and save as a new mxd. I’m not sure exactly where to start but feel there are so many layout elements that would need to be described (symbology ramp color, symbology classification, scale and extent of each layer, etc) . . . working from a template would be helpful.  The end result would be 31 different mxd and associated png.


texts: How do Natural Neighbor and IDW assign values to pixels?

texts:
According to ESRI's help files, NN interpolation calculates a value for a "query point", yet the output is ofcourse a raster.

I understand how the value of this query point is calculated based on the Voronoi diagram, but how is the eventual raster generated?

Regarding IDW, ESRI's help does mention that "The specified fixed search radius will be used for each interpolated cell (cell center) in the study area".

Should i assume that the weights, and not only the search radius itself, are calculated from the cell center?

Should i then assume that NN also refers to the centers of cells?

I have run my own experiments and my conclusion is that yes, the weights for IDW are calculated from the center of each cell, and for NN the "query point" is stationed in the center of the raster cell. Perhaps later i will add my calculations.


texts: Executing final logic after a feature is created or changed

texts:
There is a function "Python Init function" which is called when the form is opened.

Is there a function that is called when the form is closed or saved?
Or is there another way to execute some logic after I create/edited a feature?
Because I want to update the geometry of a feature based on two fields (x,y) whenever the two fields are changed.


texts:
I'm trying to run fmask on a Sentinel-2 granule.
I downloaded the 3.2 version of Windows stand alone software for Sentinel 2 from https://github.com/prs021/fmask
When I try to run fmask by calling the FmaskMSI.exe through cmd from the folder containing the Sentinel-2 Granule (or from within the granule folder) i get the following error:

No L*MTL.txt header in the current folder!
Output argument "clr_pct" (and maybe others) not assigned during call to "autoFmask".
MATLAB:unassignedOutputs

The L*MTL.txt suggests that its's looking for a Landsat type product but i specifically point to the Sentinel-2 installation of fmask

How do i get fmask to run for Sentinel-2 data?


texts: How to get a Multipolygon object from Overpass QL?

texts:
I'm new to OSM/Overpass QL and there's a simple problem that I can't solve. I've been able to figure out the a query that will return the maritime boundary of a country. For example:

(rel[boundary="administrative"][name="Puerto Rico"][admin_level="4"];>;);way(r);(._;>;);out;

This returns the boundary as a set of ways. Using overpy I've been able to get the returned ways and display them using shapely/matplotlib, and I'm confident they have the data I need.

What I can't figure out is how can I reliably convert the returned set of ways to a MultiPolygon object, ideally a shapely.geometry.MyltiPolygon.

The above query is only an example that produces a fairly simple set of ways, but there are other queries that return much more complex results, such as

(rel[boundary="administrative"]["name:en"="French Polynesia"][admin_level="3"];>;);way(r);(._;>;);out;

I understand that Overpass doesn't handle (multi)polygons, but I assume there must be a simple way to achieve this using python/GeoPandas.shapely. I've tried some manual approaches, but they are error-prone and I'd be extremely surprised if there are no implemented librariessolutions to achieve this. Any ideas?


texts: Using Make Table View tool in ArcGIS Pro?

texts:
I was going through the documentation of MakeTableView Tools here #
http://pro.arcgis.com/en/pro-app/tool-reference/data-management/make-table-view.htm

I came across the output of the tool as a out_view. The tool seems to be very helpful for my need. But when I ran that standalone from the toolbox, it successfully ran but nothing was created. I was expecting something will create. Later I realized from document that the table view that is supposed to be created by the tool is temporary and would not persist after the session ends unless the document is saved.

I am not sure what is meant by that temporary or regarding the session. If I ran the tool standalone, the session would be of its own and if I write a code then the session may be of the code to expire. Can you please help me out of this confusion Or may be you may advice me how to run the tool standalone.


texts: Using National Agriculture Statistics Service (NASS) county level crop data to estimate crop acreage in hydrologic regions?

texts:
National Agriculture Statistics Service (NASS) provides county level crop (Corn, Soybean, Sorghum, etc) data. I want to  calculate crop acres on each subwatershed from the county based data.

I'm using ArcGIS 9.3. I created a join in ArcGIS with the spreadsheet I acquired from NASS. The county level crop data is a shapefile. Watershed (subwatershed) boundaries is also in a shapefile format.

What is the best way to deal with this situation?


texts: How to display the one band layer in pseudocolor using ERDAS Imagine?

texts:
I am a new user in ERDAS imagine and i found it hard to turn the single band raster data  from grey color to pseudocolor using ERDAS imagine.

This process could be easier to be performed using ER Mapper but i can not found that there are existing color scheme available in ERDAS for better visualization purpose.


texts:
Earlier in my searches I had come across this website which provided a fast switchback between all the map providers (Apple Maps, Google, Bing, OSM) where you had the chance to compare the satellite imagery between all the different providers. Unfortunately, this website escaped my browser history and now I can't find it anywhere.

I'm aware this sounds a little delusional but is there such a site that provided this opportunity? I'm aware of geofabrik.de's comparison tool and it's similar to that but rather than side by side windows, this site was full screen.


texts: lonlat from geometry point gives strange result

texts:
I added a vector layer to my main map with the controls like the example from openlayers: http://openlayers.org/dev/examples/draw-feature.html

Now when I want to get the resulting lonlat from the point that was drawn I get strange results:

({id:"OpenLayers.Geometry.Point_81", x:6418264.3901563, y:7161906.3780079,
      bounds:{left:6418264.3901563, bottom:7161906.3780079, right:6418264.3901563,
              top:7161906.3780079}})

The main map is constructed like this:

map = new OpenLayers.Map( $('map'), {
    sphericalMercator: true,
    //maxResolution: '360/512',
    maxExtent: new OpenLayers.Bounds(-20037508.34,-20037508.34,20037508.34,20037508.34),
    projection: "EPSG:900913",
    displayProjection: new OpenLayers.Projection("EPSG:4326"),
    tileOrigin: new OpenLayers.LonLat(-180, -90),
    controls: [new OpenLayers.Control.Navigation({dragPanOptions: {enableKinetic: true}}),
                   new OpenLayers.Control.MousePosition(),new OpenLayers.Control.Zoom()],
eventListeners:{"moveend":mapEvent}
                                         });

Can somebody give me a hint, it's probably an easy thing but I'm quite new to openlayers and javascript in general :-)

Something with the conversion of projection ?


texts: How to display waterway=riverbank so that islands can be subtracted or styled separately from the water?

texts: I'd like to style my river polygons and leave the islands as background color or a different color but I can't find a way to separate the inner relations from the main riverbank multipolygon.  I've tried a bunch of k:v pairs without being able to isolate the river and I'm still new with QGIS.  What is an efficient way to do this?


texts:
When I'm doing "polygon to line" analysis on ArcGIS the result is individual line features but in QGIS is not.

How can I execute "polygon to line" in QGIS where each polygon becomes individual line feature?


texts: counting raster pixels inside or on the boundary of a postgis ingested shapefile

texts:
I'd like to count the number of raster pixels that either intersect a boundary or are contained inside the boundary and think ST_Intersects might be a good tool for this.

When using ST_Intersects on an ingested shapefile and a raster does it return true for raster pixels residing completely inside the polygon boundary or only if they overlap the border?

This is the current query I'm using but am not sure if it's picking up the pixels inside the boundary.

SELECT COUNT(*) AS intersections FROM ${rastTable} r, boundaries p WHERE p.boundaryname = '${someBoundary}' AND ST_Intersects(r.rast, p.geom);

Maybe I need to union this with an ST_Contains or is it ok as is? On further inspection I'm not even sure this is counting raster pixels. It could be counting raster tiles which isn't what I'm looking for.

EDIT

After further investigation I'm pretty sure the above query was counting raster tile intersections and not actual pixels. Here is my revised query which counts tile intersections, pixels (excluding nodatas), and pixels (including nodatas).

SELECT count(*) AS tile_intersections, SUM(ST_Count(r.rast)) as pixels_withdata, SUM(ST_Count(r.rast, false)) as total_pixels
FROM ${rastTable} r, boundaries p
WHERE p.boundaryName = ${someBoundary} AND ST_Intersects(r.rast, p.geom);

So now I just need to verify that this is including pixels contained in the boundary.

EDIT 2 (Some progress)

Messed around with this some more, the above query is also incorrect as it counts all pixels in the tiles intersecting the boundary even if they fall outside the boundary. My next attempt to solve this makes use of ST_Clip. The query does some additional work like output the result as a tiff so that I can see/debug what's happening along with using ST_SummaryStats to get other metrics I was interested in.

SELECT p.gid AS gid,
p.orgname,
ST_AsTIFF(ST_Clip(ST_Union(r.rast), p.geom, true)) AS tiffy,
ST_Count(ST_Clip(ST_Union(r.rast), p.geom, true)) AS pixel_count_excluding_nodata,
ST_Count(ST_Clip(ST_Union(r.rast), p.geom, true), false) AS pixel_count_including_nodata,
(ST_SummaryStats(ST_Clip(ST_Union(r.rast), p.geom, true), true)).*
FROM ${rastTable} r, fws_boundaries p
WHERE p.orgname = '${boundary}'
AND ST_Intersects(r.rast, p.geom)
AND r.rast_date = '${dateString}'
GROUP BY p.gid;

The above query is not optimized (I'm sure I shouldn't be clipping 4 times :) and groups by gid because some of the boundary names I use map to multiple polygons.


texts:
I have a tool  I made with a python script called 'myTool' in a toolbox called 'myBox', alias 'box'. I need to run this tool multiple times with the output of one run being one of the inputs of the next run. I would like to write a python script that would run the tool a few times. I have a script called myScript.py:

import arcpy

myTool_box('input_1.shp', 'base_case.shp', 'output_1.shp') # creates output_1.shp
myTool_box('input_2.shp', 'output_1.shp', 'output_2.shp') # creates output_2.shp
myTool_box('input_3.shp', 'output_2.shp', 'final_output.shp') # creates final_output.shp

When I run execfile(myScript.py) from the python window, I get the error:

Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "C:\Projects\myScript.py", line 3, in <module>
    arcpy.myTool_box('input_1.shp', 'base_case.shp', 'output_1.shp') # creates output_1.shp
AttributeError: 'module' object has no attribute 'myTool_box'

I assume this means that my_tool or my custom toolbox isn't on arcpy's path or I need to import the toolbox to the script, but I haven't found any information on how to tell the script where the tool resides. What am I missing here?


texts: How can I find company address with Company Name, zip codes and State of Existence?

texts:
I need to find the address of the companies located in a particular State with Company name and also zip codes.

Example:

  Company Name: Integrated Care Management

Within the following zip codes in GA -

30062, 30319, 30338, 30005, 30022, 30075, 30076, 30092, 30004

Steps:-

Currently I googling the company name and search for the address in Google maps.
Sometimes will find their company website and take the address from it. For the above example I have mentioned, I found the following address through google maps

Addresss:- 3440 Preston Ridge Road
Alpharetta, GA 30005

How can I get this done easily for more companies. Is there any website or tool to make it quite easy and less time consuming?


texts: How to compute Tree Cover Loss by year in QGIS?

texts:
I want to compute the value of tree cover loss per year from 2001 to 2014, based on the "Global Forest Change 2000 to 2014" data by Hansen et al.

I downloaded the geotif "lossyear" from http://earthenginepartners.appspot.com/science-2013-global-forest/download_v1.2.html and have added it as a raster layer to my QGIS project (version 2.16.1).

Ultimately I want to compute the tree cover loss for each year in the (green) buffer zones around household locations that I have from a survey. Setting the band rendering under properties to "Singleband pseudocolor" with equal intervals the map with the household location and buffer zones looks like this:

I used "Raster/Zonal Statistics" before to compute the % tree cover loss from a different tif that shows a dot for locations with treecover loss between 2000 and 2014 (so no distinction by year). But how can I distinguish between the different colors i.e. years in this case here, given that its just in a singleband geotif?


texts: Split lines by attributes of overlapping polygon

texts:
I have a long line, representing shoreline, with a 3km buffer overtop. The buffer has an attribute field representing the coastal characterisation for a certain length of shoreline, see image.

I am trying to split the line according to the coastal characterisation of the polygon buffer so I can get an accurate distance measure of the actual length of shoreline for each classification category


texts:
I have a shapefile of lines and another of polygons. When two or more lines intersect a polygon they should be connected (snapped), but sometimes during manual digitization this doesn't occur. I need to find and correct all these errors. I used v.clean from GRASS GIS and it automatically snapped a lot of vertices, which was good, but some of those disconnected lines are further from the tolerance I gave and if I increase the tolerance too much it begins to generate other types of errors.

Any ideas?


texts: I am new to using QGIS and for a recent project I am looking at the building heights for London. I have managed to download the free EMU Building Heights Layer and added into my existing file, through 'New Vector Layer', which displays the buildings (polygons) but struggling to get the data displayed to show different colors for the different heights.


texts: Adding raster layer in GeoServer using ArcPy gives exit code -1073741819 (0xC0000005)?

texts:
I am trying to perform a weighted overlay analysis using ArcPy. As input I'm using 5 map layers of same column and heights.

def weightOverlay():
  ambu = "rasterambu1.tif"
  fire = "rasterfire1.tif"
  hosp = "rasterhosp1.tif"
  markets = "rastermarkets1.tif"
  schools = "rasterschools1.tif"

In this way i retrieved raster layers stored in local starage(Hard disk).
Final map generated without any issues.

Then I uploaded same layers to the geoserver and tried to retrieve data as a wms.

ambu = "http://localhost:8080/geoserver/research/wms?service=WMS&version=1.1.0&request=GetMap&layers=research:ServiceAmbulance&styles=&bbox=79.8404565168211,6.862710007340316,79.89034870025498,6.980707309511519&width=330&height=768&srs=EPSG:4326&format=image%2Ftiff8"
fire = "http://localhost:8080/geoserver/research/wms?service=WMS&version=1.1.0&request=GetMap&layers=research:ServiceFireRescue&styles=&bbox=79.8404565168211,6.862710007340316,79.8903487011599,6.980707309511519&width=330&height=768&srs=EPSG:4326&format=image%2Ftiff8"
hosp = "http://localhost:8080/geoserver/research/wms?service=WMS&version=1.1.0&request=GetMap&layers=research:Hospitals&styles=&bbox=79.85332108963203,6.870961918700771,79.8883085657193,6.958916057188321&width=330&height=768&srs=EPSG:4326&format=image%2Ftiff8"
markets = "http://localhost:8080/geoserver/research/wms?service=WMS&version=1.1.0&request=GetMap&layers=research:Servicemarkets&styles=&bbox=79.8404565213477,6.862887697616292,79.8903487056314,6.980707309511519&width=330&height=768&srs=EPSG:4326&format=image%2Ftiff8"
schools = "http://localhost:8080/geoserver/research/wms?service=WMS&version=1.1.0&request=GetMap&layers=research:ServiceSchool&styles=&bbox=79.84045652315942,6.862887697616292,79.89034870025498,6.980707309511519&width=330&height=768&srs=EPSG:4326&format=image%2Ftiff8"

But this time It gives the following error

  Process finished with exit code -1073741819 (0xC0000005)

How can I get layers from GeoServer as an input for the weighted overlay analysis?


texts: Select by Location and extract features using arcpy

texts:
So I'm fairly new to arcpy and have been trying to figure out a relatively simple select by location/extract features to new layer automation in ArcMap. What I am trying to do is select features from about 40 target layers from an SDE database that are within a parish/county layer then extract them to new layers that go by the same layer name they are being outputted from but in a different location.

Below I have a couple of general issues:

1) The first for loop does not work because it will not accept "fc" as input. What I was planning out was to have it loop through each feature class and apply the parish/county layer over it to apply the selection. Of course, this does not work. This may just be due to my misuse of "fc".

2) The second for loop runs but it stops once it gets tot he second iteration/layer of the loop because it can't overwrite the filename that was already created in the first iteration. This for loop is intended to loop through each feature class and extract the same feature class to a different folder using the same layer name it was outputted from.I should also mention that the layer names go by "HM_Production.DBO.Agriculture", HM_Production.DBO.Banking", etc. and when I check the file its trying to overwrite, it has it outputted as "HM_Production.DBO". I assume this has something to do with the dot notation of the layer name?

EDIT: This is what I'm using now. It's a combination of everything so far. It goes through everything fine but is still giving me the output error of "line 2434, in CopyFeatures     raise e ExecuteError: ERROR 000210: Cannot create output C:\output\HM_Production_DBO_WATER_Water_Storage_Unit.shp Failed to execute (CopyFeatures). I should also mention that all of the scripts that have been suggested have the same issue. :

import arcpy
from arcpy import env

env.workspace = r'C:\ArcScratch\HM_Production.sde'

stateFeature = r'C:\ArcScratch\State\Parishes_LDOTD_2007.shp'
parish = "Orleans"
nameField = "PARISH"

arcpy.MakeFeatureLayer_management(stateFeature, 'SelectionStateLayer', '"' +
str(nameField) + '" = ' + "'" + str(parish) + "'")

for fc in arcpy.ListFeatureClasses():

    if arcpy.Exists('lyr'):
        arcpy.Delete_management('lyr')

    arcpy.MakeFeatureLayer_management(fc, 'lyr')
    arcpy.SelectLayerByLocation_management('lyr', "WITHIN", 'SelectionStateLayer')

    outname = fc.replace(".","_") #replace all . with _
    outpath = r"C:\output\{0}.shp".format(outname)

    if arcpy.Exists(outpath):
        arcpy.Delete_management(output)

    arcpy.CopyFeatures_management('lyr', outpath)

arcpy.Delete_management('lyr')


texts: Exporting images in native CRS and resolution

texts:
When I use “Export.image.toDrive” command in Google Earth Engine it requests to define the scale and CRS of the raster. I did a search, but could not find clear solution for it.

Is there any possibility to keep native scale and CRS of the image while exporting it?


texts: QGIS fails to clip and produces empty shapes

texts:
I've seen this type of question a few times and have tried all the different suggestions but no luck.

I have 2 layers, one polygon and one line. I want to keep the lines based on a selected polygon. I'm using the clip tool under Vector/Geoprocessing Tools/Clip. Each time, when I try to clip when setting it up and clicking OK it goes to 100% within a second or so, then the OK button stays grey for a few minutes. Eventually it comes back asking me to set the CRS for the layer but once done the layer is empty.

This is the area I'm trying to use when doing the clip. The projection for everything is NZGD2000.

I don't know if I'm doing something wrong, if the files are the problem or if there is a problem with QGIS. The files can be found at:

https://onedrive.live.com/redir?resid=85a761fa4f7e62f1!8957&authkey=!APFfcyWDxBrmy1Y&ithint=file%2czip


texts: Insert layer file to MXDs with multiple data frames using ArcPy?

texts:
I have 4  mxd's and I want to insert layer file to each mxd with position in the TOC. Each mxd has 3 and more data frame. I don't know why it does't work.
I'm trying this code:

import arcpy,os,sys
from arcpy import env

env.workspace = r"F:\Projects\Guy\mikta3\gis"
for mxdname in arcpy.ListFiles("*.mxd"):
    mxd = arcpy.mapping.MapDocument(r"F:\Projects\Guy\mikta3\gis\\" + mxdname)
    refLayer = arcpy.mapping.ListLayers(mxd, "roads", df)[0] # as display in TOC
    insertLayer = arcpy.mapping.Layer(r"F:\Projects\Guy\mikta3\gis\layers\residence.lyr")
    dfList = arcpy.mapping.ListDataFrames(mxd, "*")
    for df in dfList:
        for lyr in arcpy.mapping.ListLayers(mxd, "", df):
            arcpy.mapping.InsertLayer(df, refLayer, insertLayer, "AFTER") # BEFORE\ AFTER

    mxd.save()
del mxd


texts: Display tiles from a custom GeoTiFF overlayed on googles satellite map using leaflet

texts:
I want to overlay tiles generated from a GeoTIFF raster I created on a google satellite map. What I have tried is:

Style and output a high quality colored PNG using qgis's map composer
Tile this PNG using gdal2tiles.py
Use leaflet to point to these generated tiles

This is simple and works well, result is:

Since I want to have a "swipe" functionality to contrast my custom tiles with googles satellitle map, my tiles must be spatially aligned with googles product and have transparency where I indicate. I have not found a way to properly do this. Could someone propose a way to do this correctly?


texts: Error executing kriging in modelbuilder: 'invalid output extent'

texts:
I'm doing an interpolation using kriging in Modelbuilder, my model is this:

Every time I run the model it says:

  ERROR 010092: Invalid output extent. Failed to execute (Kriging).

Why does this happen?


texts:
I have a bunch of polygons which I loaded from a shapefile into a fields GeoDataFrame:

   geometry                                           raster_val
0  POLYGON ((8.56799 47.80800, 8.56799 47.77200, ...         1.0
1  POLYGON ((6.73199 47.12400, 6.73199 47.08800, ...         1.0

I plot them using:

fields.plot(color='white', edgecolor='black')
matplotlib.pyplot.show()

I also have a raster layer (ndvi_layer) which is basically a numpy array which I plot using:

fig, ax = matplotlib.pyplot.subplots(figsize=(10, 10))
rasterio.plot.show(ndvi_layer, ax=ax, cmap='RdYlGn')
matplotlib.pyplot.show()

Is there any way to plot the polygons on top of the raster? What if the polygons don't have the same CRS as the raster?


texts:
During the past weeks, I have tried QGIS to recognize SAGA via the processing toolbox, both being installed from the ubuntugis-unstable repository.

I have tried this on several virtual ubuntu operating systems which I specifically set up for using QGIS. However, I always get the same "Missing dependency.This algorithm cannot be run" (sic) error message. Here's some more information:

I have tried (l)ubuntu 12.04, 13.04 and 13.10 (each 64bit), installing packages from the respective repos. I
I am now an ubuntu 12.04 32bit VM
I am well aware of this post, and yes, this workaround would do the trick! However, installing from the qgis repositories is not an option for me, as I also need the orfeo toolbox packages, which are not provided there :(
Package versions are qgis 2.2.0-1, saga 2.1.1 and python-saga 2.1.1. I had the same problem with the preceding QGIS version.
I checked the PATH variable and successfully launched saga_cmd, so everything seems to be in its place.


texts:
I'm working on a webservice in ArcObjects 9.3.
It works fine when deployed on IIS, but when I try to debug some unit tests I keep getting the same error:

  The remote procedure call failed. (Exception from HRESULT: 0x800706BE)

The error occurs when calling the CreateFeatureClass(..) method throught the IDataset interface.

public void CreateOutputFeatureClass(string featureClassName)
    {
        var masterFeatureClass = Helpers.GetMaster(Session, ServerContext);
        IFeatureClassDescription fcDesc = new FeatureClassDescriptionClass();
        var ocDesc = (IObjectClassDescription)fcDesc;
        var featureDataset = featureWorkspace.OpenFeatureDataset(Session.Dataset);
        featureDataset.CreateFeatureClass(featureClassName, masterFeatureClass.Fields, ocDesc.InstanceCLSID,
            ocDesc.ClassExtensionCLSID, esriFeatureType.esriFTSimple, masterFeatureClass.ShapeFieldName, "");
    }

Based on other posts I have been reading about this error this has something to do with user permission, but I am not sure.


texts:
I want to search some keyword from GeoServer, like "E11" or "停車場" (Chinese),
My request is below:

EX: Search "E11" :

    <wfs:GetFeature service="WFS" version="1.1.0"
    outputFormat="JSON"
    xmlns:topp="http://www.openplans.org/topp"
    xmlns:wfs="http://www.opengis.net/wfs"
    xmlns:ogc="http://www.opengis.net/ogc"
    xmlns:gml="http://www.opengis.net/gml"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
    xsi:schemaLocation="http://www.opengis.net/wfs
    http://schemas.opengis.net/wfs/1.1.0/wfs.xsd">
    <wfs:Query typeName="osm_lh:planet_osm_polygon">
    <ogc:Filter>
    <ogc:PropertyIsLike wildCard="*" singleChar="#" escapeChar="!">
    <ogc:PropertyName>name</ogc:PropertyName>
    <ogc:Literal>*E11*</ogc:Literal>
    </ogc:PropertyIsLike>
    </ogc:Filter>
    </wfs:Query>
    </wfs:GetFeature>

Then I found what I want.

When I use Chinese words to search, it always return Nothing !!

EX: Search "停車場" :

    <wfs:GetFeature service="WFS" version="1.1.0"
    outputFormat="JSON"
    xmlns:topp="http://www.openplans.org/topp"
    xmlns:wfs="http://www.opengis.net/wfs"
    xmlns:ogc="http://www.opengis.net/ogc"
    xmlns:gml="http://www.opengis.net/gml"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
    xsi:schemaLocation="http://www.opengis.net/wfs
    http://schemas.opengis.net/wfs/1.1.0/wfs.xsd">
    <wfs:Query typeName="osm_lh:planet_osm_polygon">
    <ogc:Filter>
    <ogc:PropertyIsLike wildCard="*" singleChar="#" escapeChar="!">
    <ogc:PropertyName>name</ogc:PropertyName>
    <ogc:Literal>*停車場*</ogc:Literal>
    </ogc:PropertyIsLike>
    </ogc:Filter>
    </wfs:Query>
    </wfs:GetFeature>

Return:
        {"type":"FeatureCollection","totalFeatures":0,"features":[]}

I found the GeoServer log message:

    Request: getFeature
        service = WFS
        version = 1.1.0
        baseUrl = http://127.0.0.1:8090/geoserver/
        query[0]:
            filter = [ name is like *??????* ]  <==== HERE
            typeName[0] = {osm_lh}planet_osm_polygon
        outputFormat = JSON
        resultType = results
    2014-01-27 13:01:00,232 INFO [wfs.json] - about to encode JSON

So what happen to this ?? Please help !

Thanks !

Tai


texts: Using ArcPy update cursor to replace  value?

texts:
I have been combing the help associated with  (None) values in arc, and have not been able to solve the following problem:

I have a a field with mixed Null and integer values:

0

100

I would like to transform the null values to 0's in python:

   cursor = arcpy.UpdateCursor(fishJoin, ['myField']):

    for row in cursor:
        if row is None:
            row = 0
            cursor.updateRow(row)

This does not alter any of the null values.

I'm using ArcGIS Desktop 10.0.


texts: Convert KML to shapefile without losing attribute data

texts: So I want to convert KML files to SHP shapefiles to load it onto PostGIS. Currently I'm using ogr2ogr to do the conversion for me, but it wipes out all the attribute data. How do i preserve this attribute data in the resulting shapefile? Any better feature out there?


texts:
I was wondering if there is a way in OL3 or Leaflet to load another SLD file other than the one is coupled with the layer in Geoserver.
I have tried to find some topic related to this but no luck. That would be quite handy in cases you fetch data from a remote geoserver without permissions to style them.
In that case you could apply your sld style (in url format) stored locally in your server. Are you aware if this is actually supported?


texts: Choosing geospatial statistical test using building characteristics and distance together

texts:
New to geospatial analyses and I am trying to assign a building footprint to be a type of building (e.g., industrial building).

In another dataset, we have determined the average building area, average vertices, etc of industrial buildings.
Now we are trying to assign probabilities of being an industrial facility to building footprints near a GPS point (where someone said an industrial building is but it does not directly intersect with building footprint).
We assume buildings closer to the GPS point are more likely to be industrial buildings.
But we also want to match the building characteristics (area, vertices) with the known distributions. As sometimes the closest building is way different than known industrial facilities based on their area/vertices.

I've thought of building a model that predicts probability of a footprint being industrial using building characteristics, then combining with the inverse distance weight. But I'm kind of making that up and not sure how to combine.
Wondering if there are tried and true (hopefully published) methods for combining the distance from point with building/polygon characteristics?


texts:
I need GDAL with both shapefile and HDF4-Support, if possible in 64 bit.
Since the normal GDAL-versions from http://www.gisinternals.com/sdk/ have shapefile, but no HDF4 and the version from FWTools
has HDF4, but no Shapefiles (EDIT: see below), I'm afraid I have to compile a GDAL-version by myself.

The problem is that I'm 'just' a geographer with programming skills, but no real IT-expert, and I've never compiled open source by myself.

First I've downloaded the source code from http://trac.osgeo.org/gdal/wiki/DownloadSource.
The I opened makegdal90.vcproj with Visual Studio 2010. It needed to do some updates, then it opened.
If I compile it without changes, I get a GDAL-folder with a new GDAL19.dll.
If I try to use this whole folder in my C#-project that works with normal GDAL-versions, I get a error.
But if I only copy the new GDAL19.dll over the old one in the existing GDAL(32 bit)-folder it works.

Is this how it should be or is something going wrong here already ?

The second step is to switch to 64bit, in nmake.opt I changed the following entry:
 Uncomment the following if you are building for 64-bit windows
 (x64). You'll need to have PATH, INCLUDE and LIB set up for 64-bit
 compiles.
WIN64=YES
Then I changed the setting of the VS-project to X64 and compiled again.
It still works as before: I get a new GDAL19.dll that works fine in my old GDAL(64 bit)-folder.

But then I can't do the third step, including HDF4. I found and changed the following section in nmake.opt:
 Uncomment the following and update to enable NCSA HDF Release 4 support.
HDF4_PLUGIN = NO
HDF4_DIR =  C:\HDF
HDF4_LIB =  /LIBPATH:$(HDF4_DIR)\lib Ws2_32.lib
I tried it maybe 100 times, but it never worked. Sometimes I get a corrupt GDAL19.dll, sometimes I get >100 errors when compiling
and sometimes I get a working GDAL19.dll without HDF4-support. I've read almost every Google-hit in English and German I found,
but it never worked. One problem is that the files mentioned there do not exist in the HDF-versions Ive downloaded. I guess, this
tips were for much older versions of GDAL and HDF4.

So can anybody give a step by step tutorial "How to compile GDAL with HDF4 for dummies"?

I need to know which HDF4-files I need, where to get them, and what exactly I have to write in nmake.opt.

Or is there a ready to use GDAL-download with HDF4 and shapefiles somewhere ?
(32 bit would be OK to, that isn't that important)

I found out that FW-Tools-GDAL has Shapefile-support, its just my C#-project that can't use the vector-part of this old GDAL-version. In newer GDAL-versions I'm able to use vector- and rasterfiles with the same code, but unfortunately they don't have HDF4..

Maybe I need an older version of the C#-GDAL-bindings?

Has anybody used FW-Tools-GDAL in a C#-project with shapefiles?


texts: Inserting GeoJson Polygon Coordinates with MongoDB and C#?

texts:
I'm using GDAL\ogrinfo.exe to read the tab file which includes polygon/LineString/Point coordinates as a stream using the process. I'll store the coordinates in the MongoDB.

My attempt at inserting LineString was a success by using:

GeoJsonLineStringCoordinates<GeoJson2DGeographicCoordinates>

I used GeoJsonPolygonCoordinates<GeoJson2DGeographicCoordinates and kept getting the error

{"A bulk write operation resulted in one or more errors.\r\n  location object expected, location array not in correct format"}

I used the QGIS to check the tab file, if the polygon is a closed one GeoJsonPolygonCoordinates<GeoJson2DGeographicCoordinates done the job no prob. But if the tab file is like the attached image the ex-occurs.


texts: Splitting rectangles into 2 triangles using Geometry by Expression in QGIS

texts:
I am using QGIS 3.10.
I am trying to split multiple rectangle polygons diagonally into 2 parts using "Geometry by expression" under Processing.
This is for symbology only (exporting for use elsewhere), so attribute retention is not necessary.
For symbolizing diagonally split rectangles in QGIS, I successfully produced a two-layered geometry-generated symbology for each rectangle using:
Top layer:
make_triangle(point_n($geometry,1), point_n ($geometry,2),point_n( $geometry,3))
Bottom layer:
make_triangle(point_n($geometry,3), point_n ($geometry,4),point_n( $geometry,1))
This is fine for symbology, but now I need to reproduce those triangles as saved geometry.
Using Geometry by expression, I used the following, which only produces one triangle:
collect_geometries(make_triangle(point_n($geometry,1), point_n($geometry,2),point_n( $geometry,3)),(make_triangle(point_n($geometry,3), point_n ($geometry,4),point_n( $geometry,1))))
The expression builder accepts the whole expression above with no errors - just produces one element (see figure, blue polys).
Plugin Polygon divider doesn't accommodate corner to corner splits.
Also tried "Combine" expression in place of "collect geometries"
This should be simple, but I can't find any native command or plugin for the job.


texts: Seeking ArcGIS models or tools made by other people?

texts: I am new in the ArcGIS world, and I want to know if there is a site or place where people share their models and scripts for free. Like an Open source of tools for ArcGIS, for example.


texts: How does the QGIS in Boundless Desktop differ from standard QGIS?

texts:
Boundless Desktop seems to be a bundle of QGIS, PostGIS, QT Designer, a command shell, and maybe other bits. I can see that the QGIS is customised - it has some branding, and a Boundless login box. But what else is different?

The Boundless Desktop brochure is pretty light on details. My hunch (at the time of writing) is that BD is just a more attractive branding for QGIS skeptics.


texts: How to format labels with leading zeroes?

texts:
I'm trying to display a number next to a symbol in QGIS 2.0.1 - The numbers are from an integer-field in a postgresql database.

The goal is that single digit numbers get a leading 0. Examples:

12 --> 12
10 --> 10
9 --> 09
0 --> 00

In the menu

 layer properties --> labelling --> formating

I switch on "number formatting" and then enter a custom expression. I tried these two codes with various combinations of ", ' and more brackets, but none of them work:

1

CASE

      WHEN "Value" <10 THEN rpad("Value", 2, 0)

END

2

CASE

      WHEN "Value" <10 THEN "0" || "Value"

END

EDIT:

For any other person having trouble with this: Turn of "number formatting". It seems to be the reason why leading zeros are removed from the rpad function.

Enter the following in the regular "layer labeling" box:

  rpad("Field", 2, 0)

and everything will work fine. But if you also activate number formatting it will take away the leading zeros.


texts: Choropleth boundary colour change in QGIS?

texts: How do you change the colour of the boundaries between areas on a choropleth map in QGIS 2.18 Mac version?


texts:
I am new to coding and do not know what I am getting the error message: AttributeError: Object: Tool or environment  not found. My first line of code is the typical "import arcpy....." I have already tried changing my CheckOutExtension from "Spatial" to "ImageAnalyst" and got the same error message. Any suggestions?
inSegRaster = "D:\\4.Peru\\Peru_Fall\\Peru_Fall.gdb\\_2018_updatedClip"
train_features = "D:\\4.Peru\\Peru_Fall\\Peru_Fall.gdb\\just_trees"
out_definition = "D:\\4.Peru\\Peru_Fall\\SupportVector2018.ecd"
maxNumSamples = "1000"
attributes = "COLOR;MEAN"

# Check out the ArcGIS Image Analyst extension license
arcpy.CheckOutExtension("ImageAnalyst")

#Execute
arcpy.gp.TrainSupportVectorMachineClassifier(
    inSegRaster, train_features, out_definition,
    in_additional_raster, maxNumSamples, attributes)

The error message I get when running the above code is:

AttributeError: Object: Tool or environment
 not found


texts:
I am working on a large point file and am looking to get the raster value underneath each point. I am using:

UPDATE Locations a SET image3b3_sum = ST_Value(rast, 1, geom)
FROM image_sum_sum;

10 records takes 31 seconds and 1,000 takes 5:24 minutes. As I add more points
 I am noticing a linear trend on how long this will take, and for 1,000,000 points this will get nuts.

I have tried tiling the raster however it results in NULL values being assigned. What can I do to speed this up? I have looked at:
How to speed up queries for raster databases? but, as I mentioned, tiling causes null values to occur:

NOTICE:  Attempting to get pixel value with out of range raster coordinates: (-2638, -2441)
CONTEXT:  PL/pgSQL function st_value(raster,integer,geometry,boolean) line 18 at RETURN

The raster and point file have indexes. Are there constraints I can add to the raster to make it go faster?

UPDATE_1:

I have tried adding in an ST_Intersects as suggested but still yielding NULL values and long times to process:

UPDATE Locations a SET image3b3_sum = ST_Value(ST_Tile(rast,1,100,100), 1, geom)
FROM image_sum_sum b
WHERE ST_Intersects(a.geom, b.rast);


texts: Create rule-based renderer from an existing renderer using PyQGIS

texts:
I want to make a copy from an existing layer in the legend and render it with a rule-based renderer based on the original layer's graduated symbol renderer. Then I will change its rule expressions, clone all the rules (if that is even possible) and change again their expressions.

Right now I'm at step one: creating a rule-based renderer from an existing graduated symbol renderer.
With iface.actionDuplicateLayer().trigger() I clone the original layer with its renderer. Then I use convertFromRenderer function to create a QgsRuleBasedRendererV2 from the new layer's graduated symbol renderer and I assign it to the layer. But there's something I'm doing wrong because QGIS crashes at this point and I can't find the reason.

# Duplicate original layer
for layer in iface.legendInterface().layers():
    if layer.name() == 'Lines':
        iface.setActiveLayer(layer)
        iface.actionDuplicateLayer().trigger()
# Make the new layer visible and change its name
for layer in iface.legendInterface().layers():
    if layer.name() == 'Lines copy':
        vl = layer
iface.legendInterface().setLayerVisible(vl, True)
vl.setLayerName('Lines - Atlas')

# Create rule-based renderer from its graduated symbol renderer
#  and assign it to it
renderer  = vl.rendererV2()
if renderer.type() == 'graduatedSymbol':
    rSimbol = renderer.sourceSymbol()
    ruleRend = QgsRuleBasedRendererV2(rSimbol).convertFromRenderer(renderer)
    vl.setRendererV2(ruleRend) #<------- QGIS crashes here
    iface.legendInterface().refreshLayerSymbology(vl)

#Refresh canvas
canvas = iface.mapCanvas()
canvas.refresh()


texts: Show towns per province using QGIS value relation widget?

texts:
I have basic (maybe a little more) knowledge of databases and GIS.

I study Archaeology and for my thesis I am designing a local database in QGIS (consist of shapefiles with their relationships). I have to prototype the design and test it in a GIS. Since the database is for managing Dutch archaeological projectdata (not to confuse with excavation data), each shapefile contain the fieldnames "province" and "town" so we can locate the project. We can conclude that one province has many towns. I want the field "town" depend on field "province".

There are some issues... I made a .csv file with 2 fields "Value" and "Description" as like in QGIS Value Map widget. Value contains all the provinces and description the towns. In Value Relation widget in QGIS I get the province multiple times in the field. Both 'key column' and 'value column' is the value (province) of the .csv file. I would like to see each province name once.

Second, if I choose province A in field "province", I would like to see the towns of province A in field "town" and not all towns of province A,B,C...

Did I made my .csv file on a wrong way?

Or do I have to write an expression?

Which I tried, but I could not manage it to work... As I said I am an Archaeology student, not an ICT student. The first picture shows a part of my .csv in excel. The second shows what I did in the Value Relation widget. The third shows that I get all the provinces multiple times and the last shows that I get all towns (of all provinces).


texts: How to project the ecliptic on the Earth's surface in Google Earth

texts: Basically I'm trying to find all the places on Earth where the apparent path of the Sun exactly (or as nearly exact as you can get) bisects the sky or where the Sun is in the exact center of the sky at noon (relative to an observer on the Earth). This is a simple matter of obtaining the intersection of the ecliptic plane with the Earth surface, but I can't find any KMZ files that show this. Nor is there any clear method to obtain this data with Google Earth tools. I can't find any data or information on this topic at all on Google. Any leads?


texts:
How can I overlay a Vector Layer on a Raster Layer and assign NoData Values on the Raster where the Vector Layer does not cover the Raster Layer?

(preferable with ArcGIS)


texts: Which BC Albers should I use and is there any difference?

texts:
I was about to calculate the area of some polygons and before I did I did some research. I found that the area of the polygon will be represented in the units of the projection used. So I did some digging to find that the projection I need is BC Albers, but when I went into the list of projections in QGIS, I found that there were two that could the correct ones. The first is EPSG:3005 NAD83 BC ALBERS, and the second one is EPSG:3153 NAD83(CSRS)/BC ALBERS

I have projected the polygon as both and saw no difference in the result. I was looking on http://prj2epsg.org/epsg/3005 and it says:

"This CRS name may sometimes be used as an alias for NAD83(CSRS) / BC Albers. See CRS code 3153"

Just wondering if anyone knows why there are two seemingly identical projections and if there is any difference. Also what does CSRS stand for?

Thanks!


texts: Calculate list of raster from two different folder using ArcPy?

texts:
I have two different folder (workspace1 and workspace2) with Rasters data with the two different parameters. Whatever raster in workspace1 (input raster name e.g 001_mean, 002_mean, 003_mean, 004_mean ….012_mean ), I would like to execute below equation

((Raster – Min Rater value) * (0.95 – 0.01)) / (Max raster value – Min Raster Value )

and save it into outputfolder 1 name of output (e.g. 001_FPAR, 002_FPAR, 003_FPAR, 004_FPAR ….012_FPAR)

After calculating workspace1, I want to take different raster file from workspace 2 (raster name , 001_SR, 002_SR, 003_SR, 004_SR……012_SR) and multiply it with each raster file with calculated raster from workspace 1 and save it into outputfolder 2 as name of 001_APAR, 002_APAR, 003_APAR, 004_APAR ……012_APAR.

I have written the code below but its giving error.

import arcpy
arcpy.CheckOutExtension("Spatial")
from arcpy.sa import *

workspace1 = glob.glob(r"D:\MODIS-NDVI-2012\MASKED-NDVI-2012A")
workspace2 = glob.glob(r"D:\MODIS-NDVI-2012\MASKED")
outFolder1 =  r"D:\MODIS-NDVI-2012\FPAR"
outFolder2 =  r"D:\MODIS-NDVI-2012\APAR"

list = arcpy.ListRasters()

for inRaster in list:
    localRaster = Raster(inRaster)
for r in workspace1
NDVIMINResult = arcpy.GetRasterProperties_management(inRaster1, "MINIMUM")
NDVIMAXResult = arcpy.GetRasterProperties_management(inRaster1, "MAXIMUM")

NDVIMin = float(NDVIMINResult.getOutput(0))
NDVIMax = float(NDVIMAXResult.getOutput(0))
outminus = Minus(inRaster1, NDVIMin)
outminus1 = Minus(0.95-0.01)
outmultiply = Multiply(outminus,outminus1)
outminus2 = Minus(NDVIMax, NDVIMin)
outdevide = Devide(outmultiply,outminus2)
outdevide.save = (outdevide)
for r in workspace2
outmultiply1 = Multiply (outdevide,Raster)
outmultiply1.save (outmultiply1)


texts: Visualizing 3D full/solid body (like full layer) in ArcScene?

texts:
I am trying to present layer and simulated layer (how it decreased in a year, for example) in ArcScene, but showing two layers one on each other with a slight difference (like max 30 cm) does not look really great.

I tried to show the difference, but taking the heights from properties it again makes it just a layer surface, not a full body.

Is there a possibility to visualize a layer in ArcScene as a body for which you can choose the upper heights (like elevation) and the inner heights (like for ex. mean sea level or so..)?

The talk goes about Extruding, just stupidly I cannot find such option.

For illustrating. Two colorful layers (impossible to see two) represent lets say 2010 and 2011 and the difference in between them I just made a separate layer (in greyscale). And the properties box in which there is no tab 'Extrude'.


texts: Remove points which are out of shapefile (or raster) extent

texts:
My problem is quite simple. I have thousands of points (coordinates) to analyse, however some of them are so inaccurate that they "hit" the sea. I need to get rid of them. I was trying to find some R function which is capable to filter (delete) such undesirable points according to borders (mask?) of wrld_simpl (shapefile) or wrld_simpl (converted to raster)?
Please only R solution.

Here is a small working example:

# coordinates
library(maptools)
library(raster)
data(wrld_simpl)
test.points <- data.frame(lon = c(-5,-3,-4,19,19,19),
                          lat = c(45,45,45,50,51,52))

# plot (shape-file + points)
plot(wrld_simpl, axes = TRUE, col = "grey",
     xlim = c(-15,+45),
     ylim = c(+35,+60)); box()
points(test.points$lon, test.points$lat,
       col = "red", pch = 20, cex = 1)

# converting wrld_simpl to raster
ext <- extent(wrld_simpl)
xy <- abs(apply(as.matrix(bbox(ext)), 1, diff))
n <- 5
r <- raster(ext, ncol=xy[1]*n, nrow=xy[2]*n)
rr <-rasterize(wrld_simpl, r)

# plot (raster + points)
plot(rr, axes = TRUE, col = "grey",
     xlim = c(-15,+45),
     ylim = c(+35,+60)); box()
points(test.points$lon, test.points$lat,
       col = "red", pch = 20, cex = 1)


texts: Insert OLE/ActiveX objects (Excel) into ArcMap Layout using ArcObjects c#

texts:
I am trying to create an ArcMap Add-In in which I need to insert a sheet of an Excel file into the Layout using SDK.net ArcObjects. I can't find the slightest example of how to implement this. The insertion should be as via the ArcMap interface:
"Navigate to Insert > Object. In the Insert Object dialog box, select Create from file. Click Browse to go to the location of the Excel file and select it."


texts: Write numpy ndarray as TIFF with rasterio does not work

texts:
I have ndarray with shape of (12, 773, 1231) ((12 bands)).
I want to save it as array in order to reproject it. I have image with the same shape and image but reprojected so I baiscally want my ndarray "to take" the coordinates from the original image.
I am trying to use the folliwing script from this post:
img=img=rasterio.open('TOA_22_10_2020_ (1).tif')

 ##doing  processes on the image as numpy ndarray....

#now let's save it as tiff with the original image coordinates

 with rasterio.open('ff_correction.tif',
                    'w',
                    driver='GTiff',
                    height=ff_correction.shape[1],
                    width=ff_correction.shape[2],
                    count=12,
                    dtype=ff_correction.dtype,
                    crs=img.crs,
                    nodata=None, # change if data has nodata value
                    transform=img.transform) as dst:
        dst.write(ff_correction, 12)

But I recieve the next error:
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-74-786bb4d9518d> in <module>
      9                    nodata=None, # change if data has nodata value
     10                    transform=img.transform) as dst:
---> 11        dst.write(ff_correction, 1)

rasterio\_io.pyx in rasterio._io.DatasetWriterBase.write()

ValueError: Source shape (1, 12, 773, 1231) is inconsistent with given indexes 1

I don't know why the shape changes in the error.
what is the mistake?
My end goal us to reproject the ndarray to have the same coordinates as the original image, but other solutions for how to give the ndarray the original image coordinates back would work as well.
Edit:
Adding what I have done based on John answer, still getting error:
ff_correction.shape
>>>(12, 773, 1231)

np.squeeze(ff_correction).shape
>>>(12, 773, 1231)

>>> with rasterio.open('ff_correction.tif',
                    'w',
                    driver='GTiff',
                    height=ff_correction.shape[1],
                    width=ff_correction.shape[2],
                    count=12,
                    dtype=ff_correction.dtype,
                    crs=img.crs,
                    nodata=None, # change if data has nodata value
                    transform=img.transform) as dst:
        dst.write(np.squeeze(ff_correction), 12)

>>>ValueError: Source shape (1, 12, 773, 1231) is inconsistent with given indexes 1

EDIT2: When I save only one band it works. but I need all the 12...


texts: ArcMap 10.6 make column values into child group layers

texts:
ArcMap 10.6
I have a Project point feature class with the location and name of each project. I want to make a Project layer containing (sub/child) group layers for each Project name. Then I can click individual locations on and off more easily than running selections.
Right now I can manually,

Right click the data frame and make a new group layer named Projects in the data frame

For each line the the Project point feature:
2.1 Right click Project layer and create new Group Layer within, and rename this with the Project name
2.2 Create a selection of that point in the Project feature class and move it into this correctly named new Group Layer.

Is there a faster way?


texts: How change border line to no pen with python console?

texts: From Python console I load a shapefile with qgis.utils.iface.addVectorLayer. How I can remove border lines (change to no pen) in the active layer?


texts: How to insert result function in raster created pyqgis3

texts:
I am doing a script, but I have a problem with the insertion of the result that my calculation function provided in my raster output that I created. Does anyone know how to insert the output result into a raster output that has been declared.

That's what I started to do:

class ExampleProcessingAlgorithm(QgsProcessingAlgorithm):

    #Declaration de nos shapes comme des variables text
    Vector_layer_of_communes = 'Vector_layer_of_communes'
    OUTPUT_raster = 'OUTPUT_raster'
    Commune_name='Commune_name'
    INSEE_code='INSEE_code'
    EPSG_code='EPSG_code'

    def processAlgorithm(self, parameters, context, feedback):

        #Parametrage des shapes dentree
        shape_commune = self.parameterAsSource(parameters, self.Vector_layer_of_communes, context)
        #Parametrage du resultat en fichier de sorti
        output_path_raster = self.parameterAsOutputLayer(parameters, self.OUTPUT_raster, context)
        attribut_INSEE_code = self.parameterAsString(parameters, self.INSEE_code, context)
        attribut_Commune_name = self.parameterAsString(parameters, self.Commune_name, context)
        attribut_EPSG_code = self.parameterAsString(parameters, self.EPSG_code, context)

        def recuperation_cadastre_commune(layer,INSEE_code,Commune_name,EPSG_code):

            rlayer = QgsRasterLayer(urlWithParams, 'Cadastre_'+n_couche+'_'+c_insee, 'wms')
            output_path_raster=rlayer

            #My problem is here, how to insert this raster in my OUTPUT_raster created

        #Fonction de calcul du besoin
        recuperation_cadastre_commune(shape_commune,attribut_INSEE_code,attribut_Commune_name,attribut_EPSG_code)
        results = {}
        results[self.OUTPUT_raster] = output_path_raster
        return results#{}


texts: How can I create a vector layer in PyQGIS from a table without geom column?

texts:
In QGIS I've got a Spatialite database containing tables without a geometry column. With the QGIS GUI it is perfectly possible to add such a table as a 'Spatialite vector layer'.
But how to do this in Python?

On the console I've tried:

>>> uri = QgsDataSourceURI()
>>> uri.setDatabase('/<path>/test.sqlite')
>>> uri.setDataSource('', 'TestTableName','','','id')
>>> vlayer = QgsVectorLayer(uri.uri(), 'TestLayer', 'spatialite')
>>> QgsMapLayerRegistry.instance().addMapLayer(vlayer)

No result

>>> vlayer.isValid()
False

What is going wrong?

Some extra analysis:
When I add this layer by the GUI the uri-description is as follows:

>>> qgis.utils.iface.activeLayer().dataProvider().dataSourceUri()
PyQt4.QtCore.QString(u'dbname=\'/<path>/test.sqlite\' table="TestTableName" sql=')

When I add this layer by the Python console as described above, the uri-description is somewhat different:

>>> uri.uri()
PyQt4.QtCore.QString(u'dbname=\'/<path>/test.sqlite\' table="TestTableName" () sql=')

Could the empty parentheses cause the different behaviour? So maybe QGIS 1.8.0 handles empty parentheses different as no parentheses?


texts: Selecting GPS point closer to boundary of polygon in ArcGIS Desktop?

texts:
I would like to capture the gps points which are closer the polygon boundary (inside), not point in polygon.

I have arcgis arcinfo v9.3.

i can do with inner buffer of the polygon and intersection and finding the remaining gps points.  Is there any other way ?


texts: GeoServer does not render polygon properly

texts:
I am trying to a get polygon layer from CSV file. I chose to do it via rasterization and polygonization in QGIS, however, GeoServer has problems correctly rendering some polygons.

This is how the layer looks in QGIS (how it should look like):

And here is how the layer looks displayed from GeoServer (using OpenLayers):

You can see that at some point it is weirdly distorted, maybe kind of inverted. I read the Q/A GeoServer renders shape incorrectly both on Overflow and here at GIS. However, I am still lost. SRS of the layers should be set correctly as I have a lot of similar layers and its happens just for about 10% of them. If it is an issue of

  caused by how different products create polygons. (...) From the GeoServer documentation it requires that polygons have the same first and last coordinate. I would recommend that you edit the polygons in the shapefile to accomplish this.

as suggested by Darkcylde. How can I fix it? What is more important, how can I avoid it, since I am in total control of layers production and it happens for about 10 % of them?

I found out that when I rasteryze the CSV file with too high resolution causing the layers "granulation" the problem vanishes (but the picture is useless):

Here are the files which I use to create the layer. I used the modelNo of value 4.

And here is the Python code used in QGIS:

from qgis.core import *
from qgis.gui import *
from qgis.utils import *
from time import sleep
from PyQt4.QtCore import *
from PyQt4.QtGui import QImage, QColor, QPainter

# change the import
import processing

class Loader:
    def __init__(self,iface):
       """Initialize using the qgis.utils.iface
       object passed from the console.
       """
       self.iface = iface

    def create_models(self):
        print "inside a function"
        m = 11
        h = 16
        f = 4
        ss = "%02dm%02dh%dflx" % (m,h,f)

        uri = "file:///d:/GISdata/THRASH/ModGrids/MODgrid%s.dat?delimiter=%s&xField=%s&yField=%s&crs=epsg:4326" % (ss,",", "long", "lat")
        vlayer = QgsVectorLayer(uri, "MODgrid", "delimitedtext")

        writename = "d:/GISdata/THRASH/Temp/ModGridm%s.shp" % (ss)
        write = QgsVectorFileWriter.writeAsVectorFormat(vlayer,writename,"WGS84", None, "ESRI Shapefile")
        sleep(3)
        shLayer = QgsVectorLayer(writename,"PoinData", "ogr")

        raststore1 = "d:/GISdata/THRASH/Temp/testResh1.tif"
        processing.runalg("gdalogr:rasterize",shLayer,"modelNo",0,360,340,raststore1)

        NTCMname1 = "d:/GISdata/THRASH/Temp/1NTCM%s.shp" % (ss)
        processing.runalg("gdalogr:polygonize",raststore1,"DN",NTCMname1)
        player1 = QgsVectorLayer(NTCMname1, "NTCM", "ogr")
        ids = [f.id() for f in player1.getFeatures(QgsFeatureRequest().setFilterExpression(u'"DN"!=4'))]
        player1.dataProvider().deleteFeatures(ids)


texts: Color ramp that adjusts to changing scale in QGIS?

texts:
Is there a way in QGIS to have a color ramp for a raster layer (DEM layer), which adjusts to a changing scale while zooming in, and out, of the map extent?

Here is an online map made based on LIDAR, which does that: http://sggw.maps.arcgis.com/apps/webappviewer/index.html?id=c0ecc5457a2249eba0fcfaeafd40f760

When you see the map as a whole, there seems to be little variation, but once you zoom in to a larger scale, the color ramp adjusts its colors to min and max values.

It would be great to create a map with an adjusting color ramp, that way I could show the whole picture, but also export maps showing local variations.


texts: ArcObjects VB.net looping through domains of a subtype

texts:
I'm working on an ESRI-Addin written in vb.net.
I have a FeatureClass with subtypes (POI with categories). Each subtype (category) has a coded value domain defining the subcategories.
I managed to populate one combobox with the categories. See code below:
    Public Class frmFind
        Private Sub frmFind_Load(sender As Object, e As EventArgs) Handles MyBase.Load
            PopulateCategory()
        End Sub

    Private Sub PopulateCategory()
            Dim pFeatureClass As IFeatureClass = GetFeatureClass()
            Dim pSubtypes As ISubtypes = pFeatureClass
            Dim eSubtypes As IEnumSubtype = pSubtypes.Subtypes
            eSubtypes.Reset()

            Dim rCode As Integer
            Dim rName As String
            cmbCategory.Items.Clear()
            Dim dicCategories As New Dictionary(Of Integer, String)
            rName = eSubtypes.Next(rCode)
            Do While rName <> ""
                dicCategories.Add(rCode, rName)
                rName = eSubtypes.Next(rCode)
            Loop

            cmbCategory.DataSource = New Windows.Forms.BindingSource(dicCategories, Nothing)
            cmbCategory.DisplayMember = "Value"
            cmbCategory.ValueMember = "Key"
        End Sub
    ...
    End Class

Now I want to populate the second combobox with the subcategories based on the selection in the first combobox. So this is what I got so far:
    Public Class frmFind
        Private Sub frmFind_Load(sender As Object, e As EventArgs) Handles MyBase.Load
            PopulateCategory()
            PopulateSubcategory()
        End Sub

    Private Sub cmbCategory_SelectedIndexChanged(sender As Object, e As EventArgs) Handles cmbCategory.SelectedIndexChanged
        Dim pKeyValuePair As KeyValuePair(Of Integer, String) = cmbCategory.SelectedItem
        PopulateSubcategory(pKeyValuePair.Key)
    End Sub

    ...

    Private Sub PopulateSubcategory(categorycode As Integer)
        Dim pFeatureClass As IFeatureClass = GetFeatureClass()

        Dim pQFilter As IQueryFilter2 = New QueryFilter
        pQFilter.WhereClause = "KAT_ID = " & categorycode

        Dim pFeatureCursor As IFeatureCursor = pFeatureClass.Search(pQFilter, False)

        Dim pFeature As IFeature = pFeatureCursor.NextFeature

        ' Here should be a loop

        cmbSubcategory.DataSource = New Windows.Forms.BindingSource(dictSubcategory, Nothing)
        cmbSubcategory.DisplayMember = "Value"
        cmbSubcategory.ValueMember = "Key"

    End Sub

I'm struggling to get only the domains and not all features of that category. Could anybody familiar with ArcObjects give me hint?


texts: Combining text and circular background marker with ArcGIS API for Javascript?

texts:
I see that it is possible to create a circular symbol (among other shapes) with SimpleMarkerSymbol.

I also see that it is possible to create a text symbol with TextSymbol.

Is it possible to combine these two to put a bit of text, centered, over the top of a circle element?


texts: Moving points inside to polygon using PostGIS?

texts:
I have an address and building table in my PostgreSQL. There are some address points outside building polygons, the constrain is that all address must be inside of nearest polygon.

Is there a function that can help me?


texts: I need to automatically fetch Landsat images and for that it would be handy to access a ftp with images available. Does such ftp exist? Can I access it? Ultimately I want to get cloudless images.


texts: Changing the projection of the output to my OpenLayers click event handler

texts:
I've been playing around with the OpenLayers click event handler example, which you can see here.
It tell you the coordinates where you have clicked but I am trying to change this to be in decimal degrees, I believe this is epsg:900913. I've tried adding in transformation but it seems to have just made them coordinates all messed up

var tranlonlat = lonlat.transform(
                      new OpenLayers.Projection("EPSG:4326"),
                      new OpenLayers.Projection("EPSG:900913")
                      );

this is the original code:

OpenLayers.Control.Click = OpenLayers.Class(OpenLayers.Control, {
            defaultHandlerOptions: {
                'single': true,
                'double': false,
                'pixelTolerance': 0,
                'stopSingle': false,
                'stopDouble': false
            },

            initialize: function(options) {
                this.handlerOptions = OpenLayers.Util.extend(
                    {}, this.defaultHandlerOptions
                );
                OpenLayers.Control.prototype.initialize.apply(
                    this, arguments
                );
                this.handler = new OpenLayers.Handler.Click(
                    this, {
                        'click': this.trigger
                    }, this.handlerOptions
                );
            },

            trigger: function(e) {
                var lonlat = map.getLonLatFromPixel(e.xy);
                alert("You clicked near " + lonlat.lat + " N, " +
                                          + lonlat.lon + " E");
            }

        });
        var map;
        function init(){
            map = new OpenLayers.Map('map');

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

        var jpl_wms = new OpenLayers.Layer.WMS( "NASA Global Mosaic",
            "http://t1.hypercube.telascience.org/cgi-bin/landsat7",
            {layers: "landsat7"});

            jpl_wms.setVisibility(false);

            map.addLayers([ol_wms, jpl_wms]);
            map.addControl(new OpenLayers.Control.LayerSwitcher());
            // map.setCenter(new OpenLayers.LonLat(0, 0), 0);
            map.zoomToMaxExtent();

            var click = new OpenLayers.Control.Click();
            map.addControl(click);
            click.activate();

        }


texts: Filtering elements of layer with active element of QGIS atlas with geometry functions?

texts:
I'm working with QGIS 2.18

In a print composer I can filter the elements of a layer based on the attribute of the active element of an atlas. For example, show only the cities (points) of a state (polygons), when the state layer is the coverage layer of the atlas. I'm using rule-based symbolization, where both layers have a common attribute.

I would like to know if they can be filtered in the same way, but using geometry functions like the intersect. I would like to filter the elements of a layer, which intersect with the elements of the active element of the coverage layer of an atlas.

I'm trying with the intersects function, but I can not get the correct syntax.

Is it possible to make this query?


texts: What is the difference between projection and displayprojection for a Map in OpenLayers?

texts:
I am familiarizing myself with Openlayers, and am having trouble with understanding the projection and displayprojection properties for a Map.

When the map is initialized with the following code:

map = new OpenLayers.Map('map', {
        numZoomLevels: 20,
        projection: new OpenLayers.Projection("EPSG:900913"),
        displayProjection: new OpenLayers.Projection("EPSG: 4326")
        });

What exactly is happening?

The documentation of the OpenLayers.Map says that

projection -Set in the map options to override the default projection string [in] this map - also set maxExtent, maxResolution, and units if appropriate.  Default is “EPSG:4326”.

while displayProjection is the Projection used by several controls to display data to user.  If this property is set, it will be set on any control which has a null displayProjection property at the time the control is added to the map.

This does not make it clear.

Could someone explain what exactly do these two option do?


texts: How to create links between points and nearest vertex of polyline feature class?

texts:
I have a point feature class and for each point I want to find the nearest vertex on a polyline feature class. I then want to create a line between each point and it's nearest vertex on the other layer.

How is this done?

If it can't be done with a polyline layer, I can also use a point layer in lieu.


texts: Adding a mouse / touch helper to Qgis2threejs plugin

texts:
Could the Qgis2threejs plugin have a mouse / touch helper window added to the output display?

Saw this one on a Cesium display and thinking a similar helper on Qgis2threejs would make the tool more user friendly for non-GIS people.

I know working with kmz/kml files it is a pretty easy to add html linked pics for helpers, legends, title blocks, etc. and have them display on Google Earth.

You can also edit the Controls Box

This feature is available with 3DViewer(dat-gui) template.

The controls box has:
layer sub menus
Each sub menu has:

a check box to toggle layer visibility
a slider to adjust layer transparency
sub menu to control a vertically movable plane
help button to show the about box


texts: Buffering into existing FC (like in Editor menu) in Arcmap 10.3 Model Builder

texts:
I would like to buffer features into an existing FC within a GDB, similar to when buffering during an edit session by using the Buffer command in the Editor menu, not the Buffer geoprocess which creates a new FC.

Is there a way to do this in model builder?  I can only find the buffering geoprocess tool that makes a new FC.  I did consider using this, then selecting the features in the new FC and appending them into the existing FC but that adds a bunch of steps, however so far that seems like the way to get it done.

Ideally I would have a model where:

the user points to a FC to buffer
chooses the type of record in the FC to buffer
chooses which FC to buffer the features into
the model selects by attribute to focus on those records
creates buffers for the records in the existing FC
updates field values in the existing buffer FC for the new records to record what type of record they are from step 2

Any advice on that?