texts:
I started using QGIS a few weeks ago.  I've been able to create some very nice maps with it.  I've used the OpenLayers plugins to bring in different basemaps, I've brought in shapefiles and connected to our local WMS server.

But, when I connect to our local SQL server, I can bring in shapefiles, but they don't show up.  I can see all the data associated with them if I look at the properties table.  QGIS even gives me the item count, if I ask for it.

Does anyone know of a reason why I wouldn't be able to see the geometry in my maps and what I might do to fix it?

Thanks.

QGIS 2.0.1

Regarding coordinate system: The SQL database has everything in three different projections, SRID 3078, SRID 4326 and SRID 3857 (I don't know why).  I chose one set and tried each projection.  In each case, the data I brought in was the only layer in the project and the project CRS was set to match the layer CRS.  Only the data from SRID 3078 was visible.

I can't see the database directly, so I don't know if there's a geometry column.  Looking at all three datasets, each shows a feature count of 915.  I can do a query (set filter) and the same query for all three results in 116.  Nothing looks different from one to another.  It's just that two of the three don't show up in the map.

I should also note that I had our IT guy who's gathered all this data and put it in the SQL database create a shapefile from one of the datasets that I needed and couldn't read.  That showed up fine, too.


texts:
I read a lot of posts concerning that subject, nevertheless I cannot solve to problem of getting $length in the field calculator in meters rather than in decimal degress in QGIS.

The CRS of original shapefile is WGS84 (EPSG:4326). I reprojected the shape to World Euidistance Cylindrical (Sphere) (EPSG:3786). I did the same in project settings (with on the fly transformation enabled).

However I am still getting the length of the polylines in decimal degrees. Any ideas?


texts: Where is DatabaseSupport directory in the ArcGIS client?

texts:
May I know where is this folder stated in the guide?

DatabaseSupport directory in the ArcGIS client installation location and move it to a location on the Oracle server.

http://resources.arcgis.com/en/help/main/10.1/index.html#//019v0000000s000000


texts: Cannot find feature class after adding new field in ArcPy?

texts:
I'm struggling with a very basic problem. it seems though, it's a problem with arcpy itself.

I'm running the following code and it's give me a problem  in the line of:

with arcpy.da.UpdateCursor(main_points, "*" ) as cursor:

The log of the problem says is:

builtins.RuntimeError: cannot open 'main_points'

However if I remove the line of adding field:

arcpy.AddField_management(main_points, "tag_id", "SHORT", 10)

The script runs smoothly.

(PS: If I replace the line of adding field to delete field using "arcpy.DeleteField_management", the same problem happens giving me the same log error)

Here is my whole code:

import arcpy
WorkingDirectory = r"C:\Users\me\Documents"

arcpy.env.workspace = WorkingDirectory
arcpy.env.overwriteOutput=True
arcpy.CreateFileGDB_management(WorkingDirectory, "project.gdb")
arcpy.FeatureClassToGeodatabase_conversion(["main_points.shp"], "project.gdb")
arcpy.env.workspace = "project.gdb"
main_points = arcpy.ListFeatureClasses()[0]
arcpy.AddField_management(main_points, "tag_id", "SHORT", 10)
with arcpy.da.UpdateCursor(main_points, "*" ) as cursor:
    for row in cursor:
        print(row[0]) # just simple printing of Object ID
        cursor.updateRow(row)

How could I overcome this problem?


texts: How to import existing geojson data to osm and append custom tag

texts:
I have a set of geojson data where a route was recorded by a person in a wheelchair. The route goes on existing paths that were already created in OSM. How could I import this data into OSM so it will append a wheelchair accessible tag to existing paths?
Maybe another option could be that I add the tags to the geojson file and then I create my custom layer from it. The problem is that I want to use a routing service of which I have not found one that would route on a custom layer. If such a service exist please let me know.


texts: QGIS - using gdal_sieve always results in a black raster

texts:
I've got a DEM raster that I've analysed for slopes. I now need to generalise it, and it sounds like using gdal_sieve is the best way.

However, whenever I try to run it, the resulting raster is always entirely black. It seems to set each pixel to the NoValue value (in my case, -9999). My data is in floats, and goes from 0 to about 34.3 or something like that (apart from NoValues obviously).

I have used Raster>Conversion>Translate to get it to an integer, and also the Raster>Calculator to multiply the values by 100, but the results are still 0.


texts: Recommendation for a book for geographical distance calculations

texts: I am currently writing my Thesis and I need a book covering the geographical distance calculation like haversine and Vincenty. Since I am not of this field, I can't find books which might help me write. I am basically only interested in the calculations. I already have the original papers for the related stuff. But since I need some background and concept literature before I start with the calculation, I need a reference book.


texts:
I am trying to change styling of Overlay KML layer in Openlayers, but I am not able to do so.

here is my code:

<!DOCTYPE html>

<title>OpenLayers KML as Base Map Example</title>


texts: Determine the direction of groundwater flow using borehole water levels

texts:
I am working on a project that focusses on groundwater of the Marico Biosphere Reserve in South Africa and I would like to determine the general direction groundwater flows through the study area. Using the South African National Groundwater Archive, I requested information for 3000 boreholes within and surrounding the study area:

For the selected boreholes in the above image, I have the depth of the latest water level reading as well as the elevation of the borehole itself. Using either QGIS or Global Mapper, what would be the best way to determine the general flow of groundwater?


texts: How to download Launchpad ppa repositories on an android device to be installed on an airgapped ubuntu system

texts:
I have a completely air-gapped system. how can I download this(https://launchpad.net/~ubuntugis/+archive/ubuntu/ppa/+packages) repository on my android device so i can run the install command for gdal-bin

sudo apt-get install gdal-bin

need to do this because I am trying to run the command

gdalbuildvrt somegeotiff.vrt * .tif

but it says command not found, but can be installed with :

sudo apt-get install gdal-bin

to finish a vts-geospatial project.


texts:
I am storing information about a person's journey as a linestring.

Having recently discovered ST_AsSVG I tried the following query:

SELECT
ST_AsSVG(path)

FROM
    activities
WHERE
    id = 209

Which produces the path data as expexted. When i put that in to an SVG the path is exceptionally small and way off the "canvas"

See the SVG code here:

https://gist.github.com/FearMediocrity/0969b7324689a1879585

Does anyone have any suggestions to make the output more useable?


texts: Bug in ArcGIS Spatial Adjustment with Geometric Networks?

texts:
I am trying to do spatial adjustment between two datasets as shown in the the figure below. The Pipe line(selected(green color)) is part of a Geometric Network. It has three vertices.

When I perform the spatial Adjustment(tried Similarity and Affine transformation..), the Two end vertices goes where they should go but the middle vertex stays where it was before. In some other cases i observed, the entire segment stayed in its initial place without making a movement.

Is this a normal scenario while spatial adjusting features with Geometric Networks in ArcGIS? I am using ArcGIS 9.3. The spatial adjustment preview window shows data in the coorrect location as it should be.


texts:
I have a PostGIS database - using the main user I can create/update new features no worries. I have created a series of new logins and given them all full access to select/insert/update etc all the tables. However when using these new logins I get this error message when inserting new features

Could not commit changes to layer tracks

Errors: ERROR: 1 feature(s) not added.

  Provider errors:
      PostGIS error while adding features: ERROR:  permission denied for sequence TRACKS_ID_seq

The ID field is of a serial datatype using the line below to assign the value (default from PostGIS).

nextval('"TRACKS_ID_seq"'::regclass)

Not really sure what else I can try here?


texts: Copying/Pasting Information from Geoprocessing > Results > Messages

texts: I am having difficulty locating the directions for copying and pasting the information in my analysis output to a word document for future reference. Specifically, I would like to highlight and ctrl/C/V the information I attached in the screen shot.


texts: Attribute network features (lines) with data from their endpoint

texts:
Goal
To prioritize water lines in a network based on the building they service
Detail
I have a shapefile of water lines in a region. The lines are split at their vertices and form a network (more or less, will probably have to do some topo cleanup). What I want to do is find the "importance" of each pipe in the network based on the building it's connected to. That would be a simple spatial join, but I also want to attribute each connecting pipe in the network with that same info.
For example, a service line that's connected to a large building gets importance "10" (based on an attribute for that building). I want the lines that touch that service line at a shared vertex to also get a score of 10. I understand that redundant lines are present but this is for a proof-of-concept.
What I've tried
Ugly pseudocode:
for building in buildings_gdf:
    for pipe in pipes_gdf:
        if pipe.touches(building):
            pipe.buildings_connected_to = building.id + ", " # for multiple buildings

This gives me a gdf of pipes with the ID of the building they touch. But how can I find pipes that touch THESE pipes, all the way through the network, without a script that takes 10 years to run?


texts: Calculating the distance between points and multiple lines?

texts: I have a point shapefile of about 100 points and a polyline shapefile of about 200 lines and I am interested in finding the distance between every point and each individual polyline. My polyline shapefile represents individual fault lines in a region. So for each polyline (or fault line) I want to calculate the distance to each of my points. I know that I can calculate the distance between the points and the closest polyline to each point but I'm not sure how to go about calculating the distance between my points and each individual line. Any ideas?


texts: How to determine Visibility Map for Atmospheric Correction

texts:
As per Grass wiki: https://grasswiki.osgeo.org/wiki/Atmospheric_correction
the atmospheric correction could be done using a visibility raster in conjunction with an elevation raster. With regards to generating the visibility raster, I am trying to follow the stated methodology in ArcGIS (since there is no official methodology in Grass reference):
Viewshed for each pixel in raster?
(I wouldv've liked to use the Grass r.viewshed tool to stay in the same platform but it seems, it does not allow multiple coordinates of viewing position. Is creating viewshed using ArcGIS is the right approach or is there any other methodology for generating Visibility map?


texts: Insert Raster to MXD using arcpy.mapping (ArcGIS 10)

texts:
What i want to do seems simple. But again i can't even find a workaround...

I have got a lot of Raster-files (ESRI-GRID) which i need to insert into a layout mxd-file and then export the map to PDF/JPEG/whatever.

The Final  Export in Python is no Problem.
But how do I insert/replace a Raster in the mxd file?
It has to be something like:

fcs = gp.ListRasters("*", "GRID")
for fc in fcs:
    [change Raster (/ Raster Data Source?)]
    arcpy.mapping.ExportToPDF(mxd, "D:\\test.pdf")

I Hope that anybody has an idea :)
Thanks for your help!
Markus


texts: How to manipulate the lake water in the UAV image to make it better using QGIS?

texts: I am having a georeferenced UAV image over the lake. The image is a mosaic of several tiles. While doing processings by the experts the lake region ia not perfectly processed and left with voids and non-uniform appearance. How can I rectify the gap and fill the water region with uniform appearance. I have the water region digitized.


texts: How to use NOT IN correctly in a WFS CQL attribute filter?

texts:
I am using a WFS of road centrelines available here for a small mapping project in QGIS. There are a number of features in my area of interest that I just want to exclude. I figured that an attribute filter might be the easiest way to achieve this.

In the query builder dialogue within the layer properties, I added the following filter:

"t50_fid" NOT IN (3324565)

Although the "Test" button tells me that the where clause manages to return n – 1 rows, when I close the dialogue the layer fails to load, with the following exception logged:

  Layer NZ Road Centerlines (Topo, 1:50k): Download of features for
  layer data.linz.govt.nz:layer-329 failed or partially failed: Server
  generated an exception in GetFeature response:
  org.xml.sax.SAXException: Attempted to construct illegal filter - I
  dont understand the tag: fes:ValueReference. HINT: tags are
  case-sensitive! Attempted to construct illegal filter - I dont
  understand the tag: fes:ValueReference. HINT: tags are
  case-sensitive!. You may attempt reloading the layer with F5

So the NOT IN clause seems to construct an invalid CQL filter, with the layer source that is produced looking like so:

filter='"t50_fid" NOT IN (3324565)' retrictToRequestBBOX='1' srsname='EPSG:2193' typename='data.linz.govt.nz:layer-329' url='https://data.linz.govt.nz/services;key=*MYKEY*/wfs?' version='auto' table="" sql=

The following filter does work successfully, but the cost is that the syntax is harder to use if I want to exclude a longer list of features.

"t50_fid" != 3324565 (OK for one feature)

"t50_fid" != 3324565 AND "t50_fid" != 3324566 AND ... (gets cumbersome quickly, and hard to adjust).

This style of filter produces the following layer source:

filter='"t50_fid" != 3324565' retrictToRequestBBOX='1' srsname='EPSG:2193' typename='data.linz.govt.nz:layer-329' url='https://data.linz.govt.nz/services;key=*MYKEY*/wfs?' version='auto' table="" sql=

Incidentally, I figured that maybe my use of NOT IN is incorrect. "t50_fid" NOT IN (3324565,) causes no exception or warning, but fails to filter anything. "t50_fid" NOT IN [3324565] has an identical result.


texts:
I'm trying to plot a filled area from an sdf, but no matter my settings, I seem to just get an outline of the shape with a green dotted line.
I've tried with a code sample from https://developers.arcgis.com/python/guide/visualizing-data-with-the-spatially-enabled-dataframe/
from arcgis.features import FeatureLayer
fl = FeatureLayer("https://sampleserver6.arcgisonline.com/arcgis/rest/services/Census/MapServer/2")
county_sdf = fl.query("STATE_NAME='Washington'", out_sr=4326).sdf
county_sdf.head()

m9 = gis.map('Seattle, WA', zoomlevel=6)
m9.basemap='gray'
m9

county_sdf.spatial.plot(map_widget=m9,
                        cmap = 'RdPu',  # use a red to purple color map
                        symbol_type='simple',
                        symbol_style='s',
                        outline_style='s',
                        outline_color=[0,0,0,255],
                        line_width=1.0)

I get this image:

I'm looking for a filled-in shape of the state. What am I doing wrong?


texts: Select points falling outside of buffer and count

texts:
I am currently learning PostGIS and following the course offered by Datapolitan.

I want to do the steps mentioned in slide 200 https://training.datapolitan.com/qgis-training/Intermediate_GIS/#200 using a spatial query

Data can be found in https://github.com/Datapolitan-Training/qgis-training/tree/gh-pages/Intermediate_GIS/data
I used the following layers:

nyc_bike_routes_2015
nycd_16a
20150601_.........Injury.csv

The steps are pretty simple and can be done by QGIS easily as follows:

Buffer injuries by 15 feet
Select injuries near bike lanes (where the buffers intersect with the
bike lanes
Invert selection to select only injuries far from bike lane
Count off lane injuries by district

I want to do that on PostGIS and I don't seem to get the grasp of it. The query which I assume should be correct is taking too long and making my computer run out of memory

SELECT bo.*, COUNT(Q.*) as numb
FROM nyc_cd_2263 bo
LEFT OUTER JOIN (SELECT bik.boro, bik.gid, bik.instdate, inj.geom FROM injuries_2263 inj INNER JOIN bike_routes_2015 bik ON ST_DWithin(inj.geom, bik.geom, 15)) Q
ON ST_Within(Q.geom, bo.geom)
GROUP BY bo.gid

This query is intended to select and count the injuries close to the bike lanes.

Question:
I want it to do the opposite and select the injuries away from the bike lanes.

I though of using NOT ST_DWithin but this makes the query take enormous time to run.

Any idea how this could be worked out


texts: Point cluster in ArcGIS Javascript API 4.3?

texts:
I am trying to create a map with the the ArcGIS API for JavaScript v4.3. Some of the points that I want to visualize are on the same location and I would like to use some type of point clustering to render the points.

I checked the ESRI documentation for this version but I could not find something about clustering.  Is there a way to do point clustering in ArcGIS JS v4.3?


texts: How to use a property of a different interface for a created ArcObject

texts:
I have created an arc based three points (from point, through point, and to point). To create the arc I used this code:

import arcpy, comtypes.client
import comtypes.client.CreateObject as CreateObject
import comtypes.gen.esriGeometry as esriGeometry

FromObjIPoint = CreateObject(esriGeometry.Point, interface=esriGeometry.IPoint)
FromObjIPoint.PutCoords(FromX, FromY)

MidObjIPoint = CreateObject(esriGeometry.Point, interface=esriGeometry.IPoint)
MidObjIPoint.PutCoords(MidX, MidY)

ToObjIPoint = CreateObject(esriGeometry.Point, interface=esriGeometry.IPoint)
ToObjIPoint.PutCoords(ToX, ToY)

ThreePtArc = CreateObject(esriGeometry.CircularArc, interface=esriGeometry.IConstructCircularArc)
ThreePtArc.ConstructThreePoints(FromObjIPoint, MidObjIPoint, ToObjIPoint, True)

My problem is that I need to test to see the direction in which the arc is generated (clockwise or counterclockwise). There is an IsCounterClockwise property, but it is part of the ICircularArc interface, not the IConstructCircularArc interface I used to create the arc.

Any ideas on how I can test the direction of the arc?


texts: How to convert .nc4 files to .nc using R?

texts:
I have thousands of files with the .nc4 file extension and I'd like to convert them to .nc format using the programming language R.

I am using the ncdf4 package to read my data, but  I've got no clue about how I should proceed with the conversion.

Does anyone know how to do it, or where I can find helpful information on this?


texts: Keeping spatial reference using arcpy.RasterToNumPyArray?

texts:
I am using ArcGIS 10.1, and want to create a new raster based on two preexisting rasters.  The RasterToNumPyArray has a good example which I want to adapt.

import arcpy
import numpy
myArray = arcpy.RasterToNumPyArray('C:/data/inRaster')
myArraySum = myArray.sum(1)
myArraySum.shape = (myArray.shape[0],1)
myArrayPerc = (myArray * 1.0)/ myArraySum
newRaster = arcpy.NumPyArrayToRaster(myArrayPerc)
newRaster.save("C:/output/fgdb.gdb/PercentRaster")

Problem is that it strips the spatial reference and also cell size.  I figured it has to do  arcpy.env, but how do I set them based on input raster?  I cannot figure it out.

Taking Luke's answer, this is my tentative solution.

Both of Luke's solution set spatial reference, extent and cell size correctly.  But the first method did not carry data in the array correctly and output raster is filled with nodata everywhere.  His second method works mostly, but where i have big region of nodata, it fills with blocky zeros and 255s.  This may have to do with how i handled nodata cells, and i am not quite sure how i was doing it (should be another Q though).  I included images of what i am talking about.

#Setting the raster properties directly
import arcpy
import numpy

inRaster0='C:/workspace/test0.tif'
inRaster1='C:/workspace/test1.tif'
outRaster='C:/workspace/test2.tif'

dsc=arcpy.Describe(inRaster0)
sr=dsc.SpatialReference
ext=dsc.Extent
ll=arcpy.Point(ext.XMin,ext.YMin)

# sorry that i modify calculation from my original Q.
# This is what I really wanted to do, taking two uint8 rasters, calculate
# the ratio, express the results as percentage and then save it as uint8 raster.
tmp = [ np.ma.masked_greater(arcpy.RasterToNumPyArray(_), 100) for _ in inRaster0, inRaster1]
tmp = [ np.ma.masked_array(_, dtype=np.float32) for _ in tmp]
tmp = ((tmp[1] ) / tmp[0] ) * 100
tmp = np.ma.array(tmp, dtype=np.uint8)
# i actually am not sure how to properly carry the nodata back to raster...
# but that's another Q
tmp = np.ma.filled(tmp, 255)

# without this, nodata cell may be filled with zero or 255?
arcpy.env.outCoordinateSystem = sr

newRaster = arcpy.NumPyArrayToRaster(myArrayPerc,ll,dsc.meanCellWidth,dsc.meanCellHeight)

newRaster.save(outRaster)

Image showing results.  I both case nodata cells are shown yellow.

Luke's second method

My tentative method


texts:
2020-03-12T12:15:48     INFO    gdal_translate -projwin 14.780610232572156 37.735515963469275 15.197846199041757 37.43676716072688 -of GTiff cache=PreferNetwork&crs=EPSG:4326&dpiMode=7&format=GeoTIFF&identifier=1&url=http://map.sitr.regione.sicilia.it/orbs/services/raster/Slope2x2_ata0708/MapServer/WCSServer /var/folders/ds/sjyv9bzj0zg8rzpsvhbwtygr0000gn/T/processing_167546b6f1684d4a8d89e44283c0ea80/21588c8d69914a859931003cd12f1cdc/OUTPUT.tif
2020-03-12T12:15:48     INFO    GDAL execution console output
             /bin/sh: /var/folders/ds/sjyv9bzj0zg8rzpsvhbwtygr0000gn/T/processing_167546b6f1684d4a8d89e44283c0ea80/21588c8d69914a859931003cd12f1cdc/OUTPUT.tif: No such file or directory

FAILURE: No target dataset specified.
Usage: gdal_translate [--help-general] [--long-usage]
 [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/
 CInt16/CInt32/CFloat32/CFloat64}] [-strict]
 [-of format] [-b band] [-mask band] [-expand {gray|rgb|rgba}]
 [-outsize xsize[%]|0 ysize[%]|0] [-tr xres yres]
 [-r {nearest,bilinear,cubic,cubicspline,lanczos,average,mode}]
 [-unscale] [-scale[_bn] [src_min src_max [dst_min dst_max]]]* [-exponent[_bn] exp_val]*
 [-srcwin xoff yoff xsize ysize] [-epo] [-eco]
 [-projwin ulx uly lrx lry] [-projwin_srs srs_def]
 [-a_srs srs_def] [-a_ullr ulx uly lrx lry] [-a_nodata value]
 [-a_scale value] [-a_offset value]
 [-gcp pixel line easting northing [elevation]]*
 |-colorinterp{_bn} {red|green|blue|alpha|gray|undefined}]
 |-colorinterp {red|green|blue|alpha|gray|undefined},...]
 [-mo "META-TAG=VALUE"]* [-q] [-sds]
 [-co "NAME=VALUE"]* [-stats] [-norat]
 [-oo NAME=VALUE]*
 src_dataset dst_dataset


texts: gdal_translate -projwin performance slow when using VRT

texts:
I'm using gdal_translate to clip 5x5 degree tiles from a raster.

I have 4 large rasters-- one for South America, and 4 that divide Africa and Asia into 4 separate pieces. Each raster is a Byte, in WGS84, and has values of 1 (land) and NoData (ocean). I can extract a tile from South America like so:

gdal_translate -co COMPRESS=LZW --config GDAL_CACHEMAX 11% -projwin 10.0 -5.0 15.0 -10.0 1.tif my_tile.tif

This takes 40 seconds.

I also have a VRT of all the rasters together, South America and the four other tifs that compose Africa/Asia. To create the same tile, I do:

gdal_translate -co COMPRESS=LZW --config GDAL_CACHEMAX 11% -projwin 10.0 -5.0 15.0 -10.0 gadm0.vrt my_tile.tif

This takes 3 minutes.

Is this expected? Is there any way to cut down on the overhead caused by the VRT dataset? Happy to provide the actual raster data if that helps.

EDIT: See this link for the data associated data:
https://www.dropbox.com/s/y20gzv627kyjlwx/SE_Test_VRT_Data.zip?dl=0


texts: Arcgis Dynamic Mapservice from sqlserver 2008 r2

texts:
I publish the mapservices from the sqlserver2008 R2 through ArcSDE using ArcServer 10.1.
Its fine i published, now i change the attributes of my existing gdb using ArcSDE.
Now the updated information is not replicated in the mapservices.  How can i solve the issue.


texts: Creating vector joining two geopositions

texts:
I have two points in Lat, Long, Elevation format. I need to find vector/direction joining them.

Would I need to convert them both to ECEF coordinates first?


texts:
I have a geojson layer with a divIcon and now i want to change it. i have tried whit seticon but it does not work. How could i do it?

This is what I am trying:

papesjson.eachLayer(function(layer){
    if (layer.feature.properties.estado == "BAJA"){
        awesomeclase ="fa-circle-o-notch";
    }
    else if (layer.feature.properties.estado == "PENDIENTE REPOSICION"){
        awesomeclase = "fa-plus";
    }
    else if (layer.feature.properties.estado == "RETIRADA"){
        awesomeclase = "fa-times";
    }
    else {
        awesomeclase = "fa-circle";
    }
    var myIcon = L.divIcon({
        iconSize: new L.Point(50, 50),
        //iconAnchor:[50, 50],
        className: "divicon",
        html:'<i id = "'+layer.feature.properties.codigo+'" class="fa ' + awesomeclase + '" style="color:'+ getColor(layer.feature.properties.modelo) + ';">'
    });

    layer.setIcon(myIcon);
});


texts: How to display multiple plots in R Scripts running in QGIS

texts:
While plotting in R in QGIS, only the latest plot appears in 'Result viewer'. What do I do to get all the plots done during script running? Here's the code:

##showplots

bin_arr <- c("C:/Users/Madhur/Desktop/Screenshot_1.png", "C:/Users/Madhur/Desktop/Screenshot_2.png")

shp_arr <- c("C:/Users/Madhur/Desktop/SAR_Calibration/Screenshots_SAR/final.png", "C:/Users/Madhur/Desktop/SAR_Calibration/Screenshots_SAR/Screenshot_2.png")

library(imager)

im1<-load.image(bin_arr[2])
plot(im1)

im2<-load.image(shp_arr[1])
plot(im2)

im<-load.image(shp_arr[2])
plot(im)

I only get the plot corresponding to the last plot statement i.e. plot(im). I don't know if this is related or not but I ran the same script on cmd using Rscript script-name.R and not a single plot appeared.
Please help me show multiple plots in QGIS.


texts: How to remove a feature from a layer created with "Make Layer From Selected Features"?

texts:
In ArcGIS desktop, I select some features from a feature class (say with the "Select By Location"  function), and then right-click and choose "Selection--> Make Layer From Selected Features"(online help reference at bottom of page).

Now, while I am inspecting this temporary "selection layer", I see that there are some features that I really didn't intend to include and would like to remove them from this selection layer in a manual way (point-and-click).  Therefore, I would think I could select the feature (in the temporary "selection layer"), now remove it from the temporary selection layer, but cannot seem to find a way to do so.  I cannot simply start editing and delete, as this would also affect the source layer.  Is it possible to remove a feature from a selection set layer?

This is a workflow issue, as I know there are some workarounds (like export to another feature class--> select features I don't want--> and delete), but I am wondering if the direct approach is possible on the temporary selection layer.


texts: Generating flow path using r.drain from rgrass

texts:
I need to generate the flow path using r.drain from the rgrass7 package in RStudio.

I am new to GRASS GIS and prefer to use the R package instead. I want to generate a flow path from a grid cell.

General Question: For a matrix of 10*10 dimension, how to find the flow path from a given cell (not coordinate system) using GRASS GIS in R.
(https://grass.osgeo.org/grass74/manuals/r.drain.html)

Problem: Once I converted the matrix to raster, the GRASS GIS is showing an error, by stating doesn't understand S4. I do not want to use raster data; I want matrix to represent the grid cells. So that, I will generate flow path form a given grid cell - trace the lowest point from the a target cell.

rs_grass <- writeRAST(as(rs, "SpatialGridDataFrame"), "test_rs", overwrite = T)

# GRASS GIS started from R as initGRASS()

initGRASS(gisBase = "C:/Program Files/GRASS GIS 7.4.2", override = T, gisDbase = "rgrass7_TEMP")
execGRASS("g.list", parameters = list(type = "vector"))
rs_grass <- writeRAST(as(rs, "SpatialGridDataFrame"), "test_rs", overwrite = T)
execGRASS("g.region", parameters = list(vector="rs_grass", align="modis_avg@modis"))
execGRASS("v.to.rast", parameters = list(input = "rs_grass",
                                         output="rs_grass", use="cat", label_column="navn"))

Other case: Works - But has limitation for my work. I want this result from grass7

I used raster package in r - to find flow path from a grid cell.
But I want to do it in GRASS GIS in R using rgrass7 - it has more flexibility for my work.

df <- data.frame(a = c(14,5,7,9,6,14,13,12,6,8,8), b = c(9,10,10,13,7,6,7,5,9,10,11),c = c(6,8,14,15,14,12,9,5,12,14,10),
             d = c(5,8,4,13,7,5,6,9,13,10,10),e = c(13,13,10,7,9,7,10,3,13,6,11),f = c(14,15,13,7,9,13,11,14,12,12,13),g = c(6,11,9,12,13,12,5,11,8,4,6),
             h = c(15,13,14,9,7,13,14,12,4,11,11),i = c(13,15,7,8,7,11,13,8,11,7,11),
             j = c(7,9,12,15,5,7,7,5,6,15,5),k = c(14,9,6,9,9,12,15,13,14,12,7))
rs <- raster(as.matrix(df),  xmn=1, xmx=9, ymn=1, ymx=6, crs="+init=epsg:27200" )
fd <- terrain(rs, opt = "flowdir")
flowPath(fd, 49)


texts: Exporting Selected Data Driven Pages to PDF using ArcPy?

texts:
I have a map book with data driven pages that shows the fire boundaries, roads, etc for our fire rescue department. I want to automate the process of creating individual mini map books for each fire station that only include the pages within their boundary. I started off by creating a model that will iterate through the district boundary names and select all pages with those boundaries:

Now, I'd like to take the selected pages and export them to a pdf in individual folders. My goal is to have a folder for each district name and within the folder have individual pdfs of the pages that make up that district.

I'm assuming the best way to do this portion is with arcpy.mapping so I have exported the model above to a python script:

# Import arcpy module
import arcpy

# Load required toolboxes
arcpy.ImportToolbox("Model Functions")

# Local variables:
Fire_Boundaries = "Fire Boundaries_exportPDF"
Grid = "HydrantGrid_exportPDF"
Selected_District = "I_Fire Boundaries_exportPDF_NAME"

# Process: Iterate Feature Selection
arcpy.IterateFeatureSelection_mb(Fire_Boundaries, "NAME #", "false")

# Process: Select Layer By Location
arcpy.SelectLayerByLocation_management(Grid, "INTERSECT", Selected_District, "", "NEW_SELECTION")

I'm not sure where to go from here with arcpy/python?


texts: OpenLayers 3: getting rid of 'blue dot' selection icon

texts:
I am working on a digitizing tool.

There are four buttons as of now:

add point
add line
add polygon
remove geometry

The user always starts by first drawing some geometries. I am using the Draw interaction (ol.interaction.Draw) and the geometries are added to a Collection (ol.Collection).

So far so good. Now the user might decide to remove some of the drawn geometries by clicking on them.

Here is an example:

When clicked on the 'remove geometry' button, then:

drawing is removed to deactivate it
a Select interaction is created
(ol.interaction.Select) to select the geometry to remove it
when a feature is selected its ID is compared to all other features' IDs in a loop, and once the correct ID is found that feature is removed.

The logic works, but it is extremely difficult to delete a point or a line because of the blue dot that is shown when hovering, which makes it nearly impossible to click the geometries.

For example:

The following yellow point would not be removed as the blue dot is in the way:

If I move the cursor a little bit below the dot the feature will be removed:

The line is almost impossible to remove as the blue dot will move with the cursor as the cursor moves along the line:

For the polygon it works fine, as the blue dot is only shown when hovering around the edge of the geometry, but if you ignore that and click anywhere within the fill, the geometry will be removed as well:

How can I change the appearance of that blue dot? I would be OK with discarding it entirely. I have tried to add a style to the Select Interaction but that does not have any effect. It will just add another style on top of the blue dot.

EDIT: a way to solve this, perhaps, would be to look through that blue dot. Is there a way to click/select and see what is behind the blue dot?


texts: Dissolve column based on attribute in another column

texts:
I have 3 columns within a polygon layer I'm working with.

Column 1 = Landowner, 3 different values
Column 2 = Parcel, N different values
Column 3 = Management type, 5 different values

I would like to dissolve column 3 based on unique values in column 2, and again within column 1.

Example
C1 = A, B, C
C2 = 1, 2, 3, 4, 5
C3 = Z, Y, X

A     1     Z
A     1     Z
A     1     Y
A     2     X
A     2     Z
A     2     Z
B     1     Z
B     1     Y

Results
A     1     Z
A     1     Y
A     2     X
A     2     Z
B     1     Z
B     1     Y

Is there a tool in any version of ArcMap that I can set this up in or will I have to create something within the model builder or python?

Thanks


texts: Project WGS84 to Google using MapGuide and WKT

texts:
I'm having trouble projecting a simple WGS84 formatted layer over google maps. I can't find a good WKT "API" that not only shows parameters, but explains what each parameter is.

The coordinate system of the MapDefinition I am using is:

GEOGCS[
    "GCS_WGS_1984",
    DATUM[
        "D_WGS_1984",
        SPHEROID["WGS_1984",6378137.0,298.257223563]
    ],
    PRIMEM[
         "Greenwich",
         0.0
    ],
    UNIT[
        "Degree",
        0.0174532925199433
    ]
]

The PROJCS that I started writing but won't validate:

PROJCS[
    "Popular Visualisation CRS / Mercator",
    GEOGCS[
        "GCS_WGS_1984",
        DATUM[
            "D_WGS_1984",
            SPHEROID[
                "WGS_1984",
                6378137.0,
                298.257223563
            ]
        ]
        PRIMEM[
            "Greenwich",
            0.0
        ]
        UNIT[
            "Degree",
            0.0174532925199433
        ]
    ],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    AUTHORITY["EPSG","3785"],
    AXIS["X",EAST],
    AXIS["Y",NORTH]
]

For some reason mapguide won't validate this. What am I doing wrong?


texts: Syntax error: keyword can't be an expression with the function "Con"?

texts:
From How can I reclassify an integer raster to a float raster?, I would like to reclassify an integer raster "cpraster3.tif" to a float raster.
For example, I would like to reclassify "Value" = 1 into 0.27,
"Value" = 2 into 1.00 ...

Here is my code:

test_reclass = Float(Con(("cpraster3.tif", "Value" = 1, 0.27, Con("cpraster3.tif", "Value" = 2, 1.00, Con("cpraster3.tif", "Value" = 3, 2.9, Con("cpraster3.tif", "Value" = 4, 2.94, Con("cpraster3.tif", "Value" = 5, 5.37, Con("cpraster3.tif", "Value" = 6, 10.35, Con("cpraster3.tif", "Value" = 7, 16.49, Con("cpraster3.tif","Value" = 7, 22.08, 0.0)))))))))

But I have this error message:

  Failed to run script - syntax error - keyword can't be an expression


texts: Merging two neighboring cells in grid of type SpatialPolygonsDataFrame?

texts:
I have a grid of type SpatialPolygonsDataFrame (each cell has a set of data attributes). I want to merge each cell which has a surface < X with one of its neighbours. An example is depicted in the picture.

I loop over all cells and see which one has to be merged, then I get its neighbors with gTouches, but I couldn't get it work with functions like gUnaryUnion and aggregate.

Here is my code:

for(ix in(1:length(mygridshp))){
   if(mygridshp@polygons[[ix]]@area < 20000){
      neighbor_cells <- gTouches(mygridspolygons[ix], my_other_cells, byid=TRUE);
      neighbor_cells_ids <- my_other_cells$unit_id[neighbor_cells == TRUE];
      if(length(neighbor_cells_ids) > 0){
          first_neighbor <- mygridshp[mygridshp@data$unit_id==neighbor_cells_ids[1],];
          new_cell <- aggregate(rbind(first_neighbor, mygridshp[ix,]));
          new_cell <- gUnaryUnion(spgeom=union(first_neighbor,mygridshp[ix,]));
      }
   }
}

Any better idea how to do this?


texts: How to get Leaflet map.hasLayer(label) instead of map.hasLayer(layer)

texts:
I want to get label of the layer selected in option control in order to know which thematic analysis to run.
because I have just one layer, but different fillcolor regarding the label in the list control.

I want to have two separate overlays coming from the same GeoJSON source: one for population data and the other for job data? And then in the layers control select one or the other.

var mymap = L.map('mapid').setView([35.402484, 8.125763], 13);
var OpenStreetMap_Mapnik = L.tileLayer('https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
    maxZoom: 19,
    attribution: '&copy; OpenStreetMap contributors'
}).addTo(mymap);
/Color class for Population
var populationtClass = [
[0, 'rgb(249, 253, 255)'],
[500, 'rgb(195, 225, 243)'],
[1000, 'rgb(180, 219, 243)'],
[2000, 'rgb(124, 188, 227)'],
[5000, 'rgb(89, 160, 203)'],
[9000, 'rgb(30, 123, 180)'],
[19800, 'rgb(8,  86, 133)']
 ];
 //Color Class for Job
 var JobClass = [
[0, 'rgb(255, 255, 255)'],
[20, 'rgb(243, 250, 213)'],
[300, 'rgb(241, 252, 197)'],
[800, 'rgb(234, 252, 164)'],
[1000, 'rgb(224, 247, 132)'],
[2000, 'rgb(214, 250, 68)'],
[3000, 'rgb(180,222,15)']
 ];
 //function Returning the Color
 var coloring = function (thresholds, value) {
                for (var i = 0; i < thresholds.length - 1; i++) {
                    if (value < thresholds[i + 1][0])
                        return thresholds[i][1];
                }
                return thresholds[thresholds.length - 1][1];
};

//Loading Geojson of Zone1 and color it by population value or job value
var Zone1=new L.geoJson(Zone1,{
      style: function (feature) {
              var population = feature.properties.Pop_2015;
              var work = feature.properties.job_2015;
                                return {

                                    //HERE THE STUFF++++++++++++++++++++++++++++++++++++++++++++

                                    if(layer.label="Population"){
            //IF LABEL OF THE LAYER IS POP RUN POP THEMATIC COLOR
                 fillColor: coloring(populationtClass, population),
                 else if (layer.label="Jobs"){
                 //IF LABEL OF THE LAYER IS JOB RUN JOB THEMATIC COLOR
                    fillColor: coloring(JobClass, work),
                 }
                 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

                        fillOpacity: 0.75,
                        weight: 1,
                        color: 'black'
                    };
                },
onEachFeature: function (feature, layer) {
        layer.on({
                        'mousemove': function (e) {
                            e.target.setStyle({
                                weight: 7,
                                color: 'red'

                            });

                          },
                        'mouseout': function (e) {
                            Zone1.resetStyle(e.target);

                        },
                        'click': function (e) {
                            mymap.fitBounds(e.target.getBounds());

                        }});}});

L.control.layers({
'BaseMap' :OpenStreetMap_Mapnik
    },{
'Population':Zone1,
'Jobs':Zone1
}).addTo(mymap);


texts: Get and Set Lat/Long Extent of Map View in QGIS

texts:
I am new to QGIS, and assume this question is a basic one and can be easily answered. However, I was not able to find the answer by trying things out, searching the web, or reading QGIS tutorials.

I want to export a png showing a section of a map from QGIS (I used Natural_Earth_quick_start_for_QGIS.qgs with WGS 84 / Pseudo Mercator) and want to know the extent in latitude and longitude values.

I tried to do that with the Print Composer. I tried to set distinct lat/long borders in Extents, but I really wasn't able to find a mapping between lat/long values and the values in the Extents edit fields.

As stated above, I was not able to find the answer in tutorials (e.g. http://www.qgistutorials.com/en/docs/making_a_map.html) and the description in the doc isn't helpful either, see https://docs.qgis.org/2.8/en/docs/user_manual/print_composer/print_composer.html#extents

I want to be able to use two approaches:

Define a map section by mouse, as possible in the Print Composer, and obtain the corresponding lat/long borders.
Set the lat/long border values manually an obtain the corresponding map section.

How can I do that?


texts: QFileDialog.getSaveFileName() error in Plugin code for QGIS

texts:
I'm working with Plugins Builder to create a simple plugin with QGIS3 "SaveAttributes", I'm following the official tutorial on QGIS Tutorials.
I'm installed Python 3 and Python 2 is already installed with ArcGIS, I have installed also the PyQt5 and Qt Creator. (I give details because I think that the error was not provided by the code)
The error I got:

TypeError: setText(self, str): argument 1 has unexpected type 'tuple'
Traceback (most recent call last):
File "C:/Users/REDA DRISSI/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\reda\module.py", line 190, in select_output_file
self.dlg.lineEdit.setText(filename)
TypeError: setText(self, str): argument 1 has unexpected type 'tuple'

For this code:

def select_output_file(self):

    filename = QFileDialog.getSaveFileName(self.dlg, "Select output file ","", '*.txt')
    self.dlg.lineEdit.setText(filename)

I tried to use str() function :
self.dlg.lineEdit.setText(str(filename))

Gives this Error:

OSError: [Errno 22] Invalid argument: "('C:/Users/REDA DRISSI/Desktop/DATA/finaldd/reda.txt', '*.txt')"

I think the problem is related to the default version of Python used by windows and QVariant, but I don't know how to resolve the problem.


texts: What is the correct way to write the code for getting layers from the server?

texts:
I'm writing a script to call raster from the server to load them to my MXD file. It is just like when you add the data in Arcmap.

Problem is that the raster I have is on the server.

Is this correct path to retrieve them ?

raster = r"S:\Las Cruces SSO\GIS_Master_Data\NAIP_2009\Dona Ana"


texts: Data defined overrides for SVG Marker in QGIS

texts:
In QGIS 3, I have a polygon layer with numerous styles of which one is a centroid fill rendered using an SVG marker.
I would like to override the path field, so that the marker would show up only for polygons larger that x square meters.

CASE WHEN area($geometry) > 1500000
THEN 'C:/OSGEO4~1/apps/qgis/svg/gpsicons/tree.svg'
ELSE ' '
END

Yet, I am not sure what to type after ELSE, since by leaving a single blank space, I end up with question marks (unknown/undefined symbol).

Can this be done or is there a better solution?


texts:
when I customize handlerOptions of drawfeature, savestrategy does not work. Here's my code :

var handoptions = {
        handlerOptions: {
            freehand: true,
            style: {
                fillColor: "yellow",
                fillOpacity: 0.4,
                strokeColor: "green",
                strokeOpacity: 1,
                strokeWidth: 1,
                pointRadius: 6
            }
        }
    };

    draw = new OpenLayers.Control.DrawFeature(
       wfs_layer, OpenLayers.Handler.Point,handoptions,
       {
           //  title: "Draw Feature",
           //  displayClass: "olControlDrawFeaturePolygon",
           multi: true,
       }
       );

I don't know why ?


texts: How to store attribute values from a custom form into shapefile?

texts:
I have created a custom form (for mapping bird species in the field with different attributes like species, type of observation, number of birds, date, ...) with some Python logic.

One important feature of the Python logic is to provide a custom drop-down list with the possible species. A list of the most recently used species and the most often used species are added before adding the complete list of species. This way I want to minimize the amount of scrolling needed, because scrolling on a small tablet PC in the field is quite time consuming.

The form opens fine and shows the existing values (in the case of editing an existing feature). I can work with all fields as I want to, including the custom dropdown lists.

But when I finish the form by clicking the "OK" button, the changed values (in the case of editing an existing feature) or the new values (in the case of a new feature) are not written back to the attribute table of the shapefile.

So here is my question:

How can I force QGIS to write back all edited values to the attribute table??

I tried different commands like:

currentFeature.setAttributes()
currentFeature.setAttribute()
currentLayer.updateFeature(currentFeature)

but nothing has worked and I find contradictory / confusing information about these commands in the net. I guess these commands were modified in the last versions of QGIS.

System: Windows 7 64

QGIS Version: 2.8 (tried before with older versions with same result)

Python Version: 2.7


texts: Extraction climate mean values from GADM and WorldClim (Finland) using R

texts:
[Working with R 3.2.2]

I am trying to extract climate mean values country by country from GADM and WorldCLim (www.worldclim.org). The code I am using is the following one:

#Install necessary packages
install.packages('raster')
install.packages('rgdal')
install.packages('dplyr')

#Load necessary packages
library(raster)
library(rgeos)
library(dplyr)

#Set working directory
setwd('...')

#Specify target ISO country code and path to downloaded shapefile
country_name <- "FIN"
shp          <- shapefile(".../srtm_country-master/srtm/tiles.shp")       #Path to SRTM Tiles (can be found in subfolder srtm)

[Data downloaded from: https://github.com/sikli/srtm_country]

#Get country geometry first
country <- getData("GADM", country = country_name, level=2)

#Intersect country geometry with tile grid
intersects <- gIntersects(country, shp, byid=T)
tiles      <- shp[intersects[,1],]

#Download tiles
srtm_list  <- list()
for(i in 1:length(tiles)) {
  lon <- extent(tiles[i,])[1]  + (extent(tiles[i,])[2] - extent(tiles[i,])[1]) / 2
  lat <- extent(tiles[i,])[3]  + (extent(tiles[i,])[4] - extent(tiles[i,])[3]) / 2

 tile <- getData('worldclim', var='bio', res=10, lon=lon, lat=lat)

 srtm_list[[i]] <- tile
}

#Create a list with all exported .bil files from the climate data folder
ras_list <- list.files(".../wc10",full.names=TRUE, pattern=".bil$")

#Load raster mosaic into R
climatemosaic <- stack(ras_list)

#Extract mean values
assign(paste0("FIN_means"), extract(climatemosaic, country, fun=mean, na.rm=TRUE, small=TRUE, layer=1, nl=19, sp=TRUE))

#Write data.frame that keeps the needed variables
dataframe <- as.data.frame(FIN_means)

It seems to work for all the countries except FINLAND. Therefore, my guess is that the problem comes from the data from Finland. I have been trying and it does work for administrative level 0, but not for level 1 or 2 (which looks weird to me). Actually, the error message that it returns is:

Error in extent(tiles[i, ]) : error in evaluating the argument 'x' in selecting a method for function 'extent': Error in tiles[i, ] : no es posible obtener un slot ("Polygons") de un objeto de tipo "NULL"

Does anybody know if this is because the data from Finland is just available at a country level (administrative level 0)?


texts: Buffer and erase vs buffer outside_only yielding different values in ArcGIS Desktop

texts:
I am using urban atlas polygon data and wish to make a hollow buffer ring (of an arbitrary 10 m) around particular polygons.  There are two methods to do this, buffer with the setting "outside only", or buffer then erase.  I have done both of these and found that in some instances the resultant area of each method is slightly different.  I realise that buffer then erase will remove the buffer of separate polygons if they are within the buffer region, however in this example it is an isolated polygon and therefore cannot be the reason.

The figure shows an example of the polygon with the buffer.  The top attribute table (buffer_outside) is the area of the buffer by using the "outside_only" method.  The erase value is the area from using buffer then erase tools.  "Shape_Area" is an inbuilt function of the data which displays the area of the polygon and is non-editable.  To save this value, I create another column and calculate the area separately using the "calculate geometry" tool as a float.  Aside from this, all other settings are default and therefore should be equal.  Despite the funny shape of this example polygon, this issue occurs in very regular (rectangular) polygons as well.
I realise the difference is very small, but why is this happening?


texts: using QGIS 2.18.16 I am trying to figure out how to find out what polygons  in layer one (blue) intersects the polygons in layer 2 (red\pink).  I know that there will be more than one intersection for most of the polygons in layer 1 (blue) but not sure how to calc it. In a perfect world I would like to add it as a attribute to the polygon in layer 1.


texts: Which co-ordinate reference system should I use for the British Isles?

texts: I'm wondering which co-ordinate reference system would be the best to use for a project covering the British Isles, since the UK and Ireland have separate systems - British National Grid and Irish National Grid /Irish Transverse Mercator.


texts:
Shapefile has visibility attribute. This attribute is generated by algorithm and there is a manual override.

In QGIS, is there a way to apply global (per project) or local (per layer) filtering to show relevant features for given scale?

For example, if current scale is

1:25k render all features where visibility >= 5
1:10k render all features where visibility >= 6
...


texts: Join two Shapefiles based on common fields to create new polygon Shapefile?

texts:
I have two shapefiles as shown in attached picture:

Buildings footprint (polygons) which includes house number and street name beside other fields.
Red shapes (polygons) which include house number and street name with information for each building distributed in many fields.

The common fields in both shapefiles are house number and street name.

The task is creating new layer that contains Buildings footprint (target features) with corresponding Red shapes (join features) based on matching house number and street name (common fields).

For example in the attached image:

Red shape No. 1 shall be joined to building 18 (same house number & street name for both).

Red shape No. 2 shall be joined to building 16 (same house number & street name for both).

Red shape No. 3 shall be joined to building 22 (same house number & street name for both).

I tried "Spatial Join" with match options INTERSECT, HAVE_THEIR_CENTER_IN and WITHIN_A_DISTANCE, but i can't achieve the desired result.

How I would do that in ArcMap 10.3 ModelBuilder?


texts: script for automating insert scale bar, text, north arrow

texts: I've always wanted an automatic process for these tools instead of having to click insert again and again. It seems easy but I can't do it in modelbuilder and I can't seem to find anyone making a script for it on the web. I want to use a certain north arrow and scale bar each time.


texts: Simple question on my Silverlight Application.  How can I expand the Map to full extent of browser no matter the size of the browser?


texts: We can select features and display their attributes in OpenLayers. But I', curious can we add a tool similar to arcgis "Select by Attribute" in OpenLayers interface.


texts: How to add an IFeature to the Map and zoom to it ?

texts:
I have an ITable i cast it to an IFeatureLayerDefinition after casting it to an ILayer then to an IFeatureLayer, then i used my IFeatureLayerDefinition to get only one IFeature using the IFeatureLayerDefinition.DefinitionExpression so when i add my ILayer to the Map and refresh it works but i couldn't zoom to it programmatically this is my code :

private void list_commune_SelectedIndexChanged(object sender, EventArgs e)
{
    if(table_commune==null) return;

    ILayer layer = (ILayer)table_commune;
    IFeatureLayer fl = (IFeatureLayer)layer;

    object selectedCommune = ((ComboBox)sender).SelectedItem ;
    object code_commune = null;

    if (!this.commune_dic.ContainsKey(selectedCommune))
        return;
    else
        code_commune = this.commune_dic[selectedCommune];

    IFeatureLayerDefinition featureLayerDef = (IFeatureLayerDefinition)fl;
    featureLayerDef.DefinitionExpression = "NUM_COMMUNE = " + (code_commune);

    Map.AddLayer(layer);
    Map.Extent = Map.get_Layer(0).AreaOfInterest;
    Map.Refresh();
}

And this how my IFeature looks like when i added it to the Map :


texts:
I am new to qgis and trying to do 2 things. My objective is to have an accurate viewshed analysis at the end. But my DEM file is of poor resolution.

I do have a shapefile with building height details. I am trying to (Part1) find a way to make a new DEM with building heights and (Part 2) calculate the viewshed from certain buildings from the created DEM. (Part3)At a later stage, want to be able to tell from which buildings is a certain polygon on the ground visible.

Does anyone know how to do this?


texts: Extract geometries from bounding box (within and overlapping) and copy results to new empty database

texts:
I am new to PostgreSQL and GIS.

I do have a downloaded dump of OSM data (map of Germany) available, which was imported via osm2pgsql into a PostGIS-enabled PostgreSQL database germanyPG.

The query

SELECT PostGIS_Full_Version();

returns

POSTGIS="2.3.3 r15473"
GEOS="3.6.1-CAPI-1.10.1 r4317"
PROJ="Rel. 4.9.1, 04 March 2015"
GDAL="GDAL 2.2.1, released 2017/06/23"
LIBXML="2.7.8" LIBJSON="0.12"

I also have an empty PostGIS-enabled PostgreSQL database subsetGermanyPG available which has same properties (database schema, structure, tables, columns, indexes, etc.) as germanyPG database.

Question 1:

Once I've given the coordinates of a rectangle (see Figure 1).
How can I query the germanyPG database to get all objects that lie within that given rectangle or overlap with the rectangle (especially those that are much larger than the rectangle) and enclose it (see Figure 2)?

Figure 1

Figure 2

I am not sure whether ST_Intersects in combination with ST_MakeEnvelope (minLong, minLat, maxLong, max Lat, srid) will do the job because I don't have an idea how to verify the query results...

SELECT realname
FROM domain
WHERE ST_Intersects(ST_MakeEnvelope(11.5, 49.8, 12.0, 50.1, 4326),
                    geom:: GEOMETRY);

Question 2:

What is the most elegant way to copy this queried results from germanyPG database into the second database subsetGermanyPG? Notice that both databases are on the same remote host system (server side) and I need to use a JAVA JTS application (client side) who 'coordinates' respectively does the copy job.

I know that I could use JAVA to query the germanyPG database and afterward store the results of the query in a JAVA datatype (maybe one or more arrays?) and finally insert these data sequentially in subsetGermanyPG database.

But what if the query is extensive (multiple GBs) and needs a long time. I think such extreme cases would result in network traffic issues, main memory leaks, and so on.

So I hope there is a more resource-saving solution out there?

Please note, that I am only allowed to use JAVA standard library inclusive JTS.


texts: Programatically export Openlayers maps to static image file

texts:
As far as I understand this example from the OpenLayers page and the linked source code for this example, OpenLayers3 now has the offical ability to save maps as static image files.

However, the saving is done via a button on the website (and I do not fully understand how as I am not too much into Javascript).

Is there also a possibility to do this programatically in the following way?

Programatically create a webpage with an embedded OpenLayers3 map, displaying data from whatever sources
Export the webpage's canvas (all layers that would be displayed when opening the page in a browser) into a static image file

I think that creating the webpage would not be that difficult. However, the question is more on whether it is possible to use the exporting without a clickable button but by calling it via a script, command-line or similar.

I did some testing with this approach using PhantomJS to save a screenshot of the webpage. However, I have not managed yet to get this working properly (map object from my webpage not found, when found then rendering is started before all tiles have been loaded, size parameters are ignored.

Removed the OpenLayers3 dependency from the title, solutions using OpenLayers2 are also welcome.


texts: ArcGIS for Desktop exporting empty shapefile

texts:
I have a feature class in ArcGIS for Desktop 10.3. This feature class contains 12 point features.

I am attempting to export this as a shapefile;

Rightclick feature class -> Data -> Export Data (select filetype of shapefile)

However this is producing an empty shapefile. All the components of this shapefile (CPG, PRJ, SBX etc) are 1kb in size except the DBF which is 165kb and still empty (when opened in excel).

I export feature classes as shapefiles regularly and have not had this issue before.

What are the possible causes of this problem?


texts: Reporting Engines used with ArcGIS Server

texts: I'm curious what others are using to build the more elaborate reports out of ArcGIS Web API apps?  By elaborate I mean ones that include lots of dynamic charts and tables, multiple maps/data frames with legends on a single page, etc.  Seems like historically Crystal Reports and SQL Reporting Services are mentioned.  What about ReportLab?  Or with the recent improvements in ArcPy in 10.1 is it possibly better to roll something out independent of a reporting engine?


texts: I have a 2D road segment and convert it into 3D using a raster in ARCGIS 9.3. So, what i see is that the length of 3D road segment is still the same as the 2D road. How can this be? As due to elevation aspect, the 3D length should be longer than the 2D counterpart.


texts: How to Export to DXF and retain elevations of contours using QGIS?

texts: I have shape files for a whole county, and have created a base map in Qgis. What I do for my projects is zoom into the area of interest then export to dxf. The same procedure was used when I used Arcgis, however I had to go into the features and create a new feature called ELEVATION and copy the z values from the county ELEV into the newly created ELEVATION. I have tried this in Qgis but I am still not getting the contours on an elevation in cad. I have read some post about entering a code into OSGeo4WShell but the problem I have with that is my contour shape file is for a whole county and is way too big to bring into cad. I also keep getting errors so What options do I have?


texts: ArcGIS Online map labels do not appear in ArcGIS app for Android?

texts:
I have a map on ArcGIS online.  The map shows labels just fine as you can see on my web map (http://bit.ly/1q7hc0h).

The problem is that when I view the map on may Android device using the ArcGIS app the labels do not appear.

I have the most recent version of the Android app and I do not have scale dependent labels set in ArcGIS Online.

Any ideas on why the labels are not showing in the mobile version?


texts: Can an input be the same as the target with arcpy.Append?

texts:
I am basically trying to append the records of one feature class to another using arcpy.Append_management. They have similar but different schemas (they both have a uniqueid column but they also have columns that completely do not match).

osm_gdb = r"...\OSM_Differential.gdb"
osm_bing_diff_uniqueid = os.path.join(osm_gdb, "osm_bing_diff_uniqueid")
BldgFP_Parcels_20190909_gdb = "...\\BldgFP_Parcels_20190909.gdb"
osm_bing_diff_uniqueid_Layer = "osm_bing_diff_uniqueid_Layer"

arcpy.env.overwriteOutput = True
arcpy.env.workspace = BldgFP_Parcels_20190909_gdb

for fc in arcpy.ListFeatureClasses():
    state_abb = fc[:2]
    if state_abb == 'DE':
        state_feature = os.path.join(BldgFP_Parcels_20190909_gdb, fc)

        arcpy.MakeFeatureLayer_management(osm_bing_diff_uniqueid, osm_bing_diff_uniqueid_Layer, "'st_abbr' = '{}'".format(state_abb))
        arcpy.Append_management([osm_bing_diff_uniqueid_Layer, state_feature], state_feature, "NO_TEST", "", "")
        print("{} done!".format(state_abb))

Basically, I get the following error:

ExecuteError: ERROR 000572: The output cannot be the same as input.
Failed to execute (Append).

Will this just not work using Append? Should I just be using Merge instead? I am trying to keep from creating a new gdb and feature classes for the output. I really just wanted to add the rows from osm_bing_diff_uniqueid_Layer to state_feature.


texts: Length of the line inside of the polygon, ArcGIS

texts:
I need to extract (calculate) the length of the line crossing polygon. I can find the crossing points but I need 90 degrees angle from alignment projected to the edges of polygon.

Any ideas?


texts: Using processing algorithm 'native:hublines' in standalone QGIS application

texts:
I am trying to make hub lines between two vector files. I succeed manually to do that (in QGIS 3.0) and in a PyQGIS 2.18. But my code doesn't work in PyQGIS 3.0. Here it is:

output = '/home/user/Desktop'
grid = QgsVectorLayer('/home/Desktop/grid.shp',
                             "cells", "ogr")

sys.path.append('/usr/share/qgis/python/plugins')
from processing.core.Processing import Processing
Processing.initialize()
import processing

processing.algorithmHelp("qgis:hublines")
processing.run('qgis:hublines', grid, 'Province', cities, 'Province', output)

With line processing.algorithmHelp I get that hub lines doesn't exist. But it does in qgis 3.0. Also, I tried to see the list of algorithms name with processing.alglist(), but without success.

UPDATE: I tried now with 'native:hublines' like underdark said. And now I can find algorithm, but it still doesn't work and I dont have error message. Here is my new code:

output = '/home/user/Desktop/output.shp'
cities = QgsVectorLayer('path to the cities',
                             "cities", "ogr")
grid_data = QgsVectorLayer('path to the grid',"grid", "ogr")

sys.path.append('/usr/share/qgis/python/plugins')
from processing.core.Processing import Processing
Processing.initialize()
import processing

QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())

processing.algorithmHelp("native:hublines")
params = {
    'HUBS': grid_data,
    'HUB_FIELD': 'Province',
    'SPOKES': cities,
    'SPOKE_FIELD': 'Province',
    'OUTPUT': output
}

processing.run('native:hublines', params)


texts:
I am very new at this and when I georeference an image in qgis, I get the same image georeference but rotate and elongated.

I use 3 control points to georeference the image but whatever I do, I get the image with the same errors.


texts:
I was using QGIS 2.18 before the update. Adding a delimited text layer was very straightforward. However, in QGIS 3.0 I can do the same process, step by step, but it won't plot.

In summary, the link shows what I was doing in QGIS 2.18: http://www.qgistutorials.com/pt_BR/docs/importing_spreadsheets_csv.html


texts: I construct three polygons in the google map using the help polygons in google Map API for developers site. They are look like one polygon inside an another polygon .When my page is loaded all three polygon are loaded.How can i get the inner polygon when the outer polygon is zoomed or clicked in google map?


texts: estimated population around a point qgis

texts: I have a shapefile with points distributed in a polygon. The points are post office locations within a census district. How can I estimate if the post office is a representation of the concentration of population in this area recently settled (1870-1890)? Can I do any geo-statistic? Regression? Txs


texts: Change pixel values in a polygon using Google Earth Engine?

texts:
I've a classified image and polygon. The classified image pixels have two type of values (0 and 1) where 1 represents urban and 0 represents non-urban land cover. I found some wrong classified pixels as urban land and I marked the pixels' boundary by a polygon. Now I want to change the pixels' value to 0 as non-urban land cover in the polygon.

I tried to do that after converting the polygon to raster and set its value as 0. Then, I multiplied it with the main image but it didn't work.

Here is my code:

var mainImage = image1
var polygonImage = polygon.reduceToImage({
    properties: ['type'],
    reducer: ee.Reducer.first()
}).multiply(ee.Image(0));

var finalImage=mainImage.multiply(polygonImage);

I'm looking for an effective solution regarding the issue.


texts: Is there any way of making queries that contain special characters like 'ã' in spatialite-gui?

texts:
I have a Spatialite database file that has some special characters inside their columns as the following:

Here there's the word "Itapoã"... The problem is that all the times that I try to type this character in a query, spatialite-gui crashes:

The error that I see is:
An assertion failed!

../src/common/unichar.cpp(65): assert "Assert failure" failed in ToHi8bit(): character cannot be converted to single byte

I haven't found any relevant information googling this error... I know that the query works if I open the command line interface for Spatialite. The problem here is only with spatialite-gui... Is there any workaround for this kind of problem? Or I should probably just accept that inside spatialite-gui I won't be able to use special characters like 'â'?
OBS: just copying and pasting the character to the spatialite-gui interface should be enough to reproduce the problem...


texts:
I use to download layers on QGIS by Overpass using the excellent plug-in "QuickOSM".
I'm wondering why there is some wrong results, especially with points layers.
For example,  if we download the points with the key "amenity" and the value "parking" near from the Marseille's train station (France), we get some "barrier_lift_gate" which has no value for the key "amenity", and one which have the value "bicycle_parking".

Does anyone understand why ?


texts: organize the map theme and style configuration file

texts:
I am working on an android application which render some features directly, but these different features should be rendered differently at different zoom level or only be visible from a certain zoom level.

So I think I need a config file to keep this information.

During the render I will read the file, and get the matched style for a given feature, and the render it.

Since all of these are done at the client side, so I think the file size should be as less a possible(consume less memory), and the speed of querying for a matched style should be as quick as possible.

I have no experience about organize this kind of file, so I wonder if someone can give me some suggestion?

BTW,all the data no matter points or lines have an property named "uid" to identify them from others.


texts: Creating and publishing a new layer in Geoserver using Geoserver Rest API

texts:
I want to publish a new layer in Geoserver. The layer seems unpublished and I get 'Undefined' error.

const publish = function () {

return new Promise((resolve, reject) => {
    const layer = layerName
    const workspace = workspaceName
    const datastore = datastoreName

    let body = `<featureType><name>${layer}</name>
            <nativeBoundingBox>
                <minx>-180</minx>
                <maxx>180</maxx>
                <miny>-90</miny>
                <maxy>90</maxy>
                <crs>EPSG:4326</crs>
            </nativeBoundingBox>
            <latLonBoundingBox>
                <minx>-180</minx>
                <maxx>180</maxx>
                <miny>-90</miny>
                <maxy>90</maxy>
                <crs>EPSG:4326</crs>
            </latLonBoundingBox>
        </featureType>`

    console.log('body length is ', body.length);

    let options = {
        url: `http://IP:port/geoserver/rest/workspaces/datastores/featuretypes.xml`,
        method: 'POST',
        headers: {
            "Content-Type": "text.xml",
            "Content-Length": body.length
        },
        auth: {
            'user': 'user',
            'pass': 'pass'
        },

    }
    request(options, (err, resp, body) => {
        if (!err && !body) {
            return resolve()

        } else if (!err) {
            let message
            if (resp.statusCode == 400) message = "layer not found"
            else message = body

            return reject({
                status: 404,
                errors: [{
                    messages: [message]
                }]
            })

        } else {
            return reject({
                status: 500,
                errors: [{
                    messages: ["ERRO: connection with geoserver"]
                }]
            })
        }
    })
 })
}

publish().then(result => console.log(result))


texts: How do you load an asset based on the ID of the current image inside a map function?

texts:
I have a set of pre-calculated masks as image assets in a Google Earth Engine folder.  I want to map over an image collection of Sentinel images and set each image's mask based on the matching pre-calculated asset.  The precalculated mask assets are named after the Sentinel id to which they apply.

There is a similar question to this here, but my issue is that I wish to do this from inside a map function and that means that "getInfo()" will not work as you cannot use clientside operations inside a map function (see the docs).

My failing code is as follows:

var add_mask = function(image) {
  var path = ee.String("users/my_account/my_masks/");
  var mask = ee.Image(path.cat(image.id()));  //can't use getInfo() here
  return image.updateMask(mask.not());  // these masks have '1' for invalid
}

mySentinel2Collection = mySentinel2Collection.map(cloudMask);

The error I get is:

ImageCollection (Error)
Actual value for parameter 'id' must be a constant.

I have tried casting the image id as a string:

var mask = ee.Image(path.cat(ee.String(image.id())));

But this doesn't work and I get the same error.

I have also tried using combine() to join to image collections, but of course the ID don't match even though the masks use the Sentinel ID because the path to the stored image collection is part of the id.

EDIT: Actually combine() does seem to work after all - I don't know why it didn't look like it at first.  Also I don't understand why it works because the ids are demonstrably different except for the final part (after the last "/").  The documentation is opaque as to what it is actually testing here so I don't have 100% confidence - however, see my solution below.


texts: Using R and SQL together in sf - what is an efficient workflow?

texts:
Do anyone have any tips/tricks/workflows for working with an SQL query within a read_sf/st_read function call in RStudio?

(I am struggling to figure out how to make a postgrsql/postgis example reproducible)

For example I have postgresql database I am querying with something like this:

  data = read_sf(conn, query = "SELECT h.col1, h.col2, h.geom
                      FROM table AS h
                      WHERE h.col1 = \'B\'")

First, I LOVE this because it returns an sf object and I get to stay in R but leverage the database. The problem is, I am finding it messy to contruct the query in RStudio as the query string doesn't format very nicely nor does the nested quotes thing work out very well here. So my question is what are people's workflows here? Obviously query can much more complicated than this and the formatting becomes very ugly. My junky solution right now is to open a SQL editor, write the query then paste it into read_sf.

tl;dr Combining R and SQL code into read_sf presents formatting challenges that make constructing queries inefficient.


texts:
I have a simple (and maybe stupid) question for you today.

I have to compute intersection points between GIS polygons and imaginary lines I draw thanks to a script of mine... And I wonder if it would be relevant to use PostGIS instead of a Math module to catch intersection points.

Thank you.

EDIT: by math module (or library) I think about a mathematic module as I work on a Python script. So that is a "python math module".


texts: ECW not reocgnized as supported file format in GDAL operation

texts:
I've installed QGIS 3.6 for MAC, with GDAL 2.4, ECW GDAL Plugin for GDAL 2.4, and manually installed the libraries from ERDAS ECWJP2SDK 5.4.

This allows me to insert ECW files as raster layers in QGIS. I'd like to export them all as GeoTIF, so I was thinking of using GDAL Translate (convert format) in batch mode.

However, GDAL doesnt seem to recognize ECW: in QGIS messages panel I get "xyz.ecw not recognized as a supported file format.".

In QGIS preferences / GDAL I am seeing ECW with "rw+" flag and "ERDAS Compressed Wavelets (SDK 5.4)" so I think it was installed correctly (and ECW files do display, they just are not understood by GDAL).

All help deeply aprecciated.

Using macOS 10.14.3.


texts: Select_analysis Looping gives ERROR 000732

texts:
The script I use for select_analysis is no longer working. I believe part of the problem is the source data I use (and I'm not the administrator of) has changed the field properties of the attribute I use in my where_clause to string.

What I need my script to do is extract each polygon and export it to another shapefile with a naming convention that includes the where_clause field as an identifier. I use the script when the input feature has 5 or more polygons I need to extract and can be as many as 35. Also the criteria for the where_clause is numeric but not necessarily in numerical order every time. See example script below.

Right now I get the following script...

import arcpy, os

TRACT = "T407"

input_shapefile = r"S:\NRCS\HEL\Workspace\Shapefiles\T407_Flds.shp"
output_shapefile = r"S:\NRCS\HEL\Workspace\Shapefiles"

arcpy.Select_analysis('(input_shapefile)', '(output_shapefile) + (TRACT) + "_Fld5.shp"', '"clu_number" = 5')
arcpy.Select_analysis('(input_shapefile)', '(output_shapefile) + (TRACT) + "_Fld6.shp"', '"clu_number" = 6')
arcpy.Select_analysis('(input_shapefile)', '(output_shapefile) + (TRACT) + "_Fld7.shp"', '"clu_number" = 7')
arcpy.Select_analysis('(input_shapefile)', '(output_shapefile) + (TRACT) + "_Fld8.shp"', '"clu_number" = 8')
arcpy.Select_analysis('(input_shapefile)', '(output_shapefile) + (TRACT) + "_Fld16.shp"', '"clu_number" = 16')
arcpy.Select_analysis('(input_shapefile)', '(output_shapefile) + (TRACT) + "_Fld17.shp"', '"clu_number" = 17')
arcpy.Select_analysis('(input_shapefile)', '(output_shapefile) + (TRACT) + "_Fld18.shp"', '"clu_number" = 18')
arcpy.Select_analysis('(input_shapefile)', '(output_shapefile) + (TRACT) + "_Fld26.shp"', '"clu_number" = 26')
arcpy.Select_analysis('(input_shapefile)', '(output_shapefile) + (TRACT) + "_Fld38.shp"', '"clu_number" = 38')

And get the following error...

Traceback (most recent call last):
  File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\analysis.py", line 90, in Select
    raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Input Features: Dataset (input_shapefile) does not exist or is not supported
Failed to execute (Select).

So I tried using the Select tool in ArcMap and it worked extracting a single polygon, but when I copy the script and paste it to the script it doesn't like the where_clause.

Script used...

import arcpy, os

arcpy.Select_analysis(in_features="S:/NRCS/HEL/Workspace/Shapefiles/T407_Flds.shp", out_feature_class="S:/NRCS/HEL/Workspace/Shapefiles/T407_Fld5.shp", where_clause='"clu_number" = '5'')

Error Dialogue Box

There is an error in your program: invalid syntax

If I mess with the quotes around the identifier for the where_clause I get it to run but get an empty shapefile since nothing matches the clause.

Also my file paths might be part of the error, but I'm not sure. All I know about writing scripts is what I've read on the internet and it took me a long time just to get this far. Being able to loop this process would be a great time saver and why I have devoted so much time to trying to figure it out.


texts: RuntimeError: General function failure when running arcpy.da.SearchCursor

texts:
I have a script that was previously working but just started throwing an error out of the blue. Here is the traceback:

Traceback (most recent call last):

  File "<ipython-input-68-be384b1c4968>", line 4, in <module>
    nfhl_selectorz = [row[0] for row in arcpy.da.SearchCursor(S_FIRM_Pan, firm_field)]

  File "<ipython-input-68-be384b1c4968>", line 4, in <listcomp>
    nfhl_selectorz = [row[0] for row in arcpy.da.SearchCursor(S_FIRM_Pan, firm_field)]

RuntimeError: General function failure

I have checked the code at every step and the paths and workspace all seem valid. Is there anything obvious that might be causing my code to throw this error? Everything works correctly up to the line where I run SearchCursor:

import arcpy, sys, os
path = "T:/CCSI/TECH/FEMA/Datasets/NFHL/NFHL_06122018"
masterdict = {'13': ['13245'], '06': ['06007'] }

for k, v in masterdict.items():
    stateFips = k #1st iteration is '13'
    fipsList = v # first iteration is ['13245']

    for root, dirs, filename in os.walk(path):
        for dir in dirs:
            if('NFHL_' + stateFips in dir and '.gdb' in dir):
                temp_path = os.path.join(root, dir)
                arcpy.env.workspace = temp_path

                featureclasses = arcpy.ListFeatureClasses()
                for fc in featureclasses:
                    if fc == 'S_FIRM_Pan':
                        S_FIRM_Pan = fc
                        arcpy.env.workspace = S_FIRM_Pan
                        field_names = [f.name for f in arcpy.ListFields(S_FIRM_Pan)]
                        for fip in fipsList:
                            fipC = fip + 'C'
                            firm_field = "DFIRM_ID"
                            nfhl_selectorz = [row[0] for row in arcpy.da.SearchCursor(S_FIRM_Pan, firm_field)]


texts: Is there a QGIS plugin/inbuilt tool that allows the user to export a georeferenced image into a SQL blob field?


texts: How to Insert Feature into Feature Class (in Memory) ESRI

texts:
I'm looking for suggestions on how to insert a new feature instance into a feature class that is in memory -- assuming that I don't have access to the feature workspace. For example, the feature class has 4 "rows", and I want to add another to make the total 5.

So something like this ....

void InsertFeatureInstance(IFeatureClass fc, IFeature feature)
{
     // Do the work here.
}

I have a good idea how to do this when I have the feature workspace, using a feature cursor, calling CreateFeature, etc. But in this case, I want to be able to accomplish it all without touching the underlying geodb.

Working in a stand-alone C# application, ArcGIS 10, and ArcGIS Engine.

--

EDIT: Here is the solution I came up with. Some of it was copied from other code I got from somewhere else (can't remember where).

public static void InsertFeatureInstance(IFeatureClass targetFeatureClass, IFeature sourceFeature)
{
    IFeature newfeature = targetFeatureClass.CreateFeature();
    var simplifiedFeature = newfeature as IFeatureSimplify;
    IGeometry myGeometry = sourceFeature.ShapeCopy;
    simplifiedFeature.SimplifyGeometry(myGeometry);
    newfeature.Shape = myGeometry;

    for (int i = 1; i < sourceFeature.Fields.FieldCount; i++)
    {
        IFeatureClass sourceFeatureClass = (IFeatureClass) sourceFeature.Class;

        string fieldName = sourceFeature.Fields.Field[i].Name;
        bool bCondition1 = fieldName == sourceFeatureClass.ShapeFieldName;
        bool bCondition2 = (sourceFeatureClass.LengthField != null &&
                            fieldName == sourceFeatureClass.LengthField.Name);
        bool bCondition3 = (sourceFeatureClass.AreaField != null &&
                            fieldName == sourceFeatureClass.AreaField.Name);

        if (!(bCondition1 || bCondition2 || bCondition3))
        // Don't do shape fields
        {
            int myTargetFieldId = targetFeatureClass.FindField(fieldName);
            // Id of field in source feature
            newfeature.Value[myTargetFieldId] = sourceFeature.Value[i]; // Copy value
        }
    }

    newfeature.Store();
}


texts: How to bulk print maps using ArcGIS or QGIS?

texts: I have a polygon layer that I would like to print images or PDFs from in bulk as there are a large number of them. Rather than having to manually scroll and print each individual'site' in the print composer (QGIS 2.0 or ArcGIS 10) are there any tools that I could use or scripts etc?


texts: Postgis how to add a vertex to a linestring at intersection with a point?

texts:
I'm using postgis 2.0 and I want to add a vertex to a linestring (pipes) at intersection position with a point (startpoint of waterlateral).
Can anyone help?


texts:
I want to select entries of my GPS-Point data only recorded from 7 a.m. to 7 p.m. and close out all others. Right now, its in a String format (00:00:06).
I tried "time": >= '07:00:00' and "time": <= '19:00:00' ", but it doesnt work. I think, its because of the wrong format of the column. Maybe it works with an integer or a daytime format? I tried to convert, but i did not find a way until now.

I am using Qgis 3.0.3.


texts: Wireframe display of terrain in ArcGlobe?

texts:
I have successfully imported a terrain dataset into ArcGlobe which provides elevation values to the globe surface.

Is there an easy way I can now change the symbology/display of the terrain/globe surface to a 'wireframe' style display, ideally with customisable grid spacing etc? Or would I need to drape a grid image over the terrain surface?


texts: Getting data closest to specified point in a grib file (Python)?

texts:
I have a grib file. I would like to be able to do something like this:

grib_data = grib_data_nearest_point(grib_file='grib_file.grib2', lat=15.5, lon=23.3)

Is there a library that can accomplish this in python, or some library that has a workaround that can accomplish this? I am currently reading in every point in the grib file, and performing the distance calculation on every single point, which is incredibly slow.

Here is my comparison of my old method and the new method suggested:

import timeit
import pygrib
import math
import gdal
import numpy as np
from geopy.distance import geodesic

class SpeedGribOld:

    def get_closest_datapoint(self, location_to_match, grib_file):

        grib_file_handler = pygrib.open(grib_file)
        data = grib_file_handler[1].data()[0].flatten()
        grib_lats = grib_file_handler[1].data()[1].flatten()
        grib_lons = grib_file_handler[1].data()[2].flatten()

        match_lat = location_to_match[0]
        match_lon = location_to_match[1]

        distances = []
        grib_indices = []
        for index, grib_lat in enumerate(grib_lats):
            grib_lon = grib_lons[index]

            # Only check grib points within 2 degrees lat/lon of the powerplant
            if abs(grib_lat - match_lat) < 2:
                if abs(grib_lon - match_lon) < 2:
                    distances.append(geodesic((grib_lat, grib_lon), (match_lat, match_lon)))
                    grib_indices.append(index)
        # Get min distance index
        index_min_dist = np.argmin(distances)
        closest_data = data[grib_indices[index_min_dist]]
        print("OLD: Closest data: {}".format(closest_data))

class SpeedGribNew:
    def get_closest_datapoint(self, location_to_match, grib_file, band_i=1):
        """
        Gets the correspondent value to a latitude-longitude pair of coordinates in
        a grib file.

        :param grib_file:   path to the grib file in disk
        :param match_lat:         latitude
        :param match_lon:         longitude
        :param band_i:      band index
        :return:            scalar
        """

        match_lat = location_to_match[0]
        match_lon = location_to_match[1]

        # open the grib file, get the specified band and read it as an array
        ds = gdal.Open(grib_file, 0)
        band = ds.GetRasterBand(band_i)
        arr = band.ReadAsArray()

        # get origin's coordinates, pixel width and pixel height
        # the GetGeoTransform method returns the skew in the x and y axis but you
        # can ignore these values
        ox, pw, xskew, oy, yskew, ph = ds.GetGeoTransform()

        # calculate the indices (row and column)
        i = math.floor((oy - match_lat) / ph)
        j = math.floor((match_lon - ox) / pw)

        # close the file
        del ds, band

        # index the array to return the correspondent value
        closest_data = arr[i, j]
        print("NEW: Closest data: {}".format(closest_data))

grib_files = ['/path/to/GFS_UREL.grib']

start = timeit.default_timer()
SpeedGribOld().get_closest_datapoint(
    location_to_match = np.array([10, 30]),
    grib_file = grib_files[0]
)
stop = timeit.default_timer()
print('Old Time: ', stop - start)

start = timeit.default_timer()
SpeedGribNew().get_closest_datapoint(
    location_to_match = np.array([10, 30]),
    grib_file = grib_files[0]
)
stop = timeit.default_timer()
print('New Time: ', stop - start)

Output:

OLD: Closest data: 2.2722020149230957
Old Time:  2.053238352003973
NEW: Closest data: -3.877797842025757
New Time:  0.03821680700639263


texts: Shapefile Conversion to ERDAS Signature editor

texts: I'm trying to run a supervised classification for an image. I was working with it in ArcMap and created some training data. I wanted to see if I could get a better result with Erdas Imagine using the same training data. I was able to convert the original training data from ArcMap to an AOI in Erdas, but can't seem to go from there to the signature editor so I can run the supervised classification. I'd like to use the same training data in Erdas that I used in ArcMap- is that possible or should I just start the classification from scratch?


texts:
If I try to load the following data using the Open Oracle Layer in QGIS, the layer is disabled and QGIS won't load it :

insert into rLINES( id, geom) select 1, SDO_GEOMETRY( 2002, NULL, NULL, SDO_ELEM_INFO_ARRAY( 1,2,1),
        SDO_ORDINATE_ARRAY( 195421.075730808, 326086.472701877,
                            195431.899003386, 326086.909972049,
                            195441.359009782, 326088.651032466 )) from dual;

If I use this data, QGIS will display it :

insert into RLINES( id, geom) select 1, SDO_GEOMETRY( 2002, NULL, NULL, SDO_ELEM_INFO_ARRAY( 1,2,1),
        SDO_ORDINATE_ARRAY( 10, 10,
                            15, 15,
                            20, 20 )) from dual;

Why is that ?

If I exit QGis and then start it again, the data is displayed.


texts: r - create multipolygon from overlapping polygons using sf package

texts:
I'm trying to remove overlapping areas from multiple polygons using sf like the following:

# sample polygon
poly <- data.frame(
    lon = c(0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 2, 2, 0.8, 1, 1, 2, 2, 1, 1),
    lat = c(0, 0, 1, 1.5, 0, 1, 1, 2, 2, 1, 1, 1, 2, 2, 1, 0, 0, 1, 1, 0),
    var = c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3 ,3 ,3 ,3 ,3, 4 ,4 ,4, 4, 4)
  ) %>%
  st_as_sf(coords = c("lon", "lat"), dim = "XY") %>% st_set_crs(4326) %>%
group_by(var) %>%
  summarise(geometry = st_union(geometry), do_union = F) %>% st_cast("POLYGON")

In my data there are several hundred polygons and many of them are overlapping multiple times. So it would be great if there is a way to remove overlapping areas from a sf dataframe like in the example. I tried something with st_difference but I thought there might be a better way than applying this to each pair of polygons which would result in too many permutations.


texts: Using touch/tap features in mobile devices with OpenLayers?

texts:
I'm new on OpenLayers and I'm developing an app for my final project career.

I develop a map viewer where you can tap the features displayed on it and a popup is shown.
When I try the app in PC browser using click you can tap the feature with no difficulties.

The problem is when I want to tap one of this features in the mobile version. It's very difficult to tap and hit the feature.

My code:

        map.on('click'/*'touch'*/, function(evt) {

            popup.hide();
            popup.setOffset([0, 0]);

            var feature = map.forEachFeatureAtPixel(evt.pixel, function(feature) {

                if (feature.getGeometry().getType() == 'Point'){
                return feature;}

            });

         });

I tried with different possible solutions but any of them solved the issue  100%:

var snap = new ol.interaction.Snap({
            source: source.getSource(),
            pixelTolerance: 30
            });
            map.addInteraction(snap);

or

hitTolerance: 15

or

getclosestfeature();

Do you know how to resolve it?


texts: Running PyQGIS standalone script on Ubuntu terminal (docker)

texts:
I am running a PyQGIS standalone script on Ubuntu terminal in a QGIS docker:

from qgis.core import *
QgsApplication.setPrefixPath("/path/to/qgis/installation", True)
qgs = QgsApplication([], False)

But I get the following error:

  QStandardPaths: XDG_RUNTIME_DIR not set, defaulting to
  '/tmp/runtime-root' qt.qpa.screen: QXcbConnection: Could not connect
  to display Could not connect to any X display.


texts:
The public organization I work for has some archived hyperspectral data (.raw - radiance, non georectified).
Has anyone processed hyperspectral data for consumption within ArcMap? If so, what sort of products have you created from the hyperspectral data?
I would like to see this information utilized somehow instead of sitting on the server collecting dust, so to speak.


texts: How to display multiple records in single feature point using ArcMap 10.x?

texts:
I have a list of touristic sights in a point layer (around 500 records). Each record has a spatial location field - municipality (town, city). In many cases there are more than one sights in a single municipality. I need to display them in a map of 1 : 100 000 scale.
So, logically, if I displayed them in a normal way, the dots would be overlapping (one city - many dots).

Is there any feasible and cartographically common way to do this?
What came to my mind:

having only one dot in the map which would say upon clicking/hovering on it there are many records in that spot
maybe there's a setting that will keep the dots from overlapping, perhaps spreading them in a circle or something along this way?

(I'm not sure how to correctly state the question in the title, please, if you can edit it so it is more understandable, do so)


texts: Inserting value in middle of attribute table in ArcGIS Desktop?

texts:
I using this attribute table of a point layer:

and i want to insert a new row between lines 7 and 8 when i start editing the layer, and not as The default- the new feature will be added automatically to the end of the attribute table. My goal is to join this attribute table with another attribute table from another layer.
How can i do it in ArcGIS Desktop 10.3?


texts: Shortest distance for multiple starting points

texts:
I am looking to calculate the shortest distance along a road network for multiple starting points. The reason I am not doing this manually is that I have more than 5000 starting points.
Does anyone know how to do this in QGIS 3? I have researched extensively with no results for this particular problem.

Input data:

Snapshot of the results, which are probably wrong:


texts:
I started with a linear model (lm) and I want to make some spatials models so I tried first a Moran test to lm.

Can we reject the null hypothesis with a p value of 0.11 in order to continue working and trying spatial models?

I have tried several specifications for lms and all of them have bigger p values.

I want just to justify why I will try spatial models.


texts: "Not a line layer:" error message when using mmqgis animate lines

texts:
I'm working on animating lines that grow to its  full length and this error message

  Not a line layer:

keeps popping up while attempting to use the MMQGIS function.

File used is a feature line generated from ArcMap. MMQGIS plugin was run in QGIS 2.8.3.

Just wondering if I've missed some vital step(s) or do I have to convert my feature line (shp file) into another format.


texts: Reproject file to line up with geographic features on a map

texts:
I have been trying to reproject data to match the "geographic features" on a map (and also shapefiles) in QGIS or on the web. I'm ingesting raw GRIB data (available here). I'm using python to open the grib and get the values. I'm the using gdal to try to save to file.
Here is the python script:
import numpy as np
import pygrib
from osgeo import osr, gdal

grbs = pygrib.open('data/hrrr.t23z.wrfprsf00.grib2')
lats, lons = grbs[607].latlons()
surfacePres = grbs[607].values

flipped_data = np.flipud(surfacePres)

nRows = surfacePres.shape[0]
nCols = surfacePres.shape[1]

xres = (lats.max() - lats.min()) / nRows
yres = (lons.max() - lons.min()) / nCols

driver = gdal.GetDriverByName("GTiff")
dst_ds = driver.Create('sfcpres.tif', nCols, nRows, 1, gdal.GDT_Float32)
dst_ds.SetGeoTransform([lons.min(), yres, 0, lats.max(), 0, -xres])

srs = osr.SpatialReference()
srs.SetWellKnownGeogCS("WGS84")
dst_ds.SetProjection(srs.ExportToWkt())

band = dst_ds.GetRasterBand(1)
band.SetNoDataValue(9999)
band.WriteArray(flipped_data)

I have printed nRows, nCols, xres, and yres.
nRows 1059
nCols 1799
xres 0.029723824652320985
yres 0.04067720231374391

When I go to visualize my tif on a map (opening in QGIS), the file is not matching up with the base maps. You can clearly see this near the great lakes, and if you look closely on the west coast, the "darker colors" (purple, etc) are over water, when they should actually be over land.

When I go to open the original grib2 file (attached above - the "file"), QGIS uses proj4 in order to overlay it on my map. This depiction is what I am expecting.

Finally, this last image shows the difference between the two. The gray image is what I am currently generating, and the colored image is what I am trying to get.


texts:
I am working with open data from the City of Toronto. I am using a base layer of property boundaries which can be found here:

http://www1.toronto.ca/wps/portal/contentonly?vgnextoid=c3a11a23871b9310VgnVCM1000003dd60f89RCRD&vgnextchannel=1a66e03bb8d1e310VgnVCM10000071d60f89RCRD

The default projection system I am using is WGS84.  However, when I place the map, at any scale, in print composer the scale bar is completely wrong.  I would like to find a way to measure distance from lot lines accurately, but need an accurate scale.  I have tried changing the CRS settings without luck.

Also, is there a way in QGIS like there is in Autocad where if I know the distance between two points and input that, the scale of the shapefile will adjust to my measurement accordingly?

I am using QGIS 2.2 btw.


texts: When I use the ogr2ogr -f command to create a geopackage the first layer imports correctly but every layer after that drops fields. Can anybody explain why it does this?


texts:
I am trying to write some mocha/expect unit tests for my OpenLayers based app. I tried using jsdom to simulate the DOM, and ended up having to add also xmldom, canvas, and now it still won't work since it requires  requestAnimationFrame. So maybe I should use a different approach.

I am now looking at using phantomjs. Are there any examples for testing ol3 in that way? is it feasible? Or maybe I should just test it in a normal browser?


texts: ArcPad 10.0.2 - Removing shortcuts from ArcPad today

texts: I have created a shortcut that goes to ArcPad today menu, but now I want to remove it.  How does one do this?


texts: Defining Max Value for Breaks in Symbology instead of data driven Breaks

texts:
I'm mapping positivity rates per Zip code. I have four main race counts in my attribute table and I duplicate this shapefile in order to create 4 different map layouts; White positivity per Zip code, Black positivity per Zip code, Hispanic positivity per Zip code and Asian positivity per Zip code.
When I use graduated symbology and classify my values by selecting the related field on my attribute table, QGIS sets my max value to that field's max value. So when I put these into my layout it looks like each race has the maximum positivity in some areas while in reality their positivity rate is much lower.
I'm trying to map these values (rates) in a way where the max value is 100 across all 4 layers representing %100 positivity and map races accordingly.
I.e if in zip code 55555 whites have %20 blacks have %30 hispanics have %40 and asians have %50 and when I classify this layers in whites' layer my color that represent the maximum value (let's say red) shows up completely red but the same red is assigned to %30 for blacks, %40 for hispanics and %40 for asians. But I'm trying to map these so %20 for whites, blacks, hispanics and asians are exactly the same color.
Is there a way to achieve this without changing classification numbers by hand every time I click on classify?