texts: What is the proper "title" property for the Terrain with Labels Basemap?

texts:
EDIT: Actually this posted link ended up helping me quite a bit but don't mark it as duplicate because the problem I had using it stemmed from needing to define the map before defining the basemapGallery. See ADDENDUM below.

Note that this post does NOT help me ...
How do I use the Terrain with Labels basemap in the Map constructor

I am using PHP and trying to pass in this variable (among others) representing the basemap id "basemap_5" to determine the selected basemap when a new page is "launched" (basically my page launches a new map with the same center (XY) and zoom level).

Then using $_GET I can access the php variable in my javascript to determine which basemap was being viewed and use those variables like this ...

map = new Map("map", {
basemap: bMapName,
center: passedCenter,
zoom: zoomLevel

});

I added a new argument for the basemap (as shown above ... "bMapName") and they all work except this one (basemap_5 aka Terrain with Labels). The error I get is ...

Map.setBasemap: Unable to find basemap definition for: "Terrain with Labels". Try one of these: streets,satellite,hybrid,topo,gray,oceans,national-geographic,osm

I sure wish the error message were updated to include the terrain with labels map!

I also tried "World Terrain Base" since that appears to be the name of the service
Link to Service

EDIT:
Let me add to this that I'm using the BasemapGallery widget. The problem centers around needing to identify the basemap before creating the basemapGallery (which has "map" set for it's map property). If I can just get around that, I think I could get this to work. I have tried the samples below and from the post I linked to above. I can add custom maps but they aren't recognized (obviously) before they have been added. So how can I get around that?

    map = new Map("map", {
    basemap: bMapName,
    center: passedCenter,
    zoom: zoomLevel
  });

  //add the basemap gallery, in this case we'll display maps from ArcGIS.com including bing maps
  basemapGallery = new BasemapGallery({
    showArcGISBasemaps: true,
    map: map
  }, "basemapGallery");

ADDENDUM
I got this to work by defining "map" conditionally and setting the basemap property to a created basemap (as in the posted link) in the condition that one of the two maps that weren't working (my custom basemap and the "Terrain with Labels" basemap) were the ones being referenced "over the wire". setting the basemap property as shown in the linked SE post was actually the most helpful part of this but the answer provided here was similarly helpful.

  switch (bMap) {
  case "basemap_0":
    bMapName = "MD Imagery";
    break;
  case "basemap_1":
    bMapName = "osm";
    break;
  // SEVERAL MORE CASES HERE; TRUNCATED FOR SE POST ...
  }
mdImagelayer = new ArcGISTiledMapServiceLayer("http://geodata.md.gov/imap/rest/services/Imagery/MD_SixInchImagery/MapServer");
  if (bMapName === "terrain") {
    map = new Map("map", {
      basemap: new Basemap({id: 'terrain',
                  layers: [new BasemapLayer({url: 'http://services.arcgisonline.com/arcgis/rest/services/World_Terrain_Base/MapServer'}),
                    new BasemapLayer({url: 'http://services.arcgisonline.com/ArcGIS/rest/services/Reference/World_Reference_Overlay/MapServer'})
                  ]
            }),
      center: passedCenter, // [-79.2, 39.5]
      zoom: zoomLevel // 12
    });
  } else if (bMapName === "MD Imagery") {
    map = new Map("map", {
      basemap: new Basemap({id: 'MD Imagery', layers: [mdImagelayer]}),
      center: passedCenter, // [-79.2, 39.5]
      zoom: zoomLevel // 12
    });
  } else {
  // You may wish to change the id to map or mapDiv (if that is the map you are using
    map = new Map("map", {
      basemap: bMapName,
      center: passedCenter, // [-79.2, 39.5]
      zoom: zoomLevel // 12
    });
  };

  //add the basemap gallery, in this case we'll display maps from ArcGIS.com including bing maps
  basemapGallery = new BasemapGallery({
    showArcGISBasemaps: true,
    map: map
  }, "basemapGallery");


texts: How to extract the US from a North America Shapefile?

texts:
I have a north america shapefile, downloaded from
ftp://ftp.epa.gov/wed/ecoregions/cec_na/NA_CEC_Eco_Level2.zip

How can I get a new shapefile from this showing only us? Is there a way doing this in QGIS?

Thanks


texts:
I am trying to do kriging with my data. But there is a problem here. I want to use krige command in gstat package. But I do not have "Spatial object with prediction/simulation locations".

How I can make SpatialPolygons?

I have a shapefile but the problem is that I want to do kriging outside of my shapefile.


texts: GRASS GIS: creating raster with a certain gradient and extent

texts:
Constraints
- Use only GRASS GIS
- Automated using pyGRASS

Input

DEM
River raster OR river as vector (both are computed)
gradient of the valley (as a float)
vector line with direction and magnitude of the gradient

The idea
I want to create a script (pyGRASS) to execute the following steps automatically.

The general idea is to construct a "virtual surface" layer. This surface layer will be a raster with a certain extent and gradient. The direction of this gradient will be described by a vector line.

Thus, in essence I want to create a raster (of certain width) along a vector line with a certain gradient. Thus the raster will start with a certain height which will decrease along the line up untill the end of the line.

Creating "virtual surface"

I am not sure how I will actually create this surface. I've looked into r.gradient but here I have to specify the direction of the gradient by hand. I already have the direction and magnitude (vector)

Has anyone an idea or a different approach to solve this problem? Faced it before?


texts:
I have a very tricky problem:
I have a Postgis 2.0 database, used as datastore for Geoserver 2.1.3

I have a table with a geometry column (named shape), with SRID 4326

I have created a view like this:

CREATE VIEW pippo AS SELECT st_geometryn(shape,1) as geom FROM events

And then I published the view using geoserver. Everything works, except for GetFeatureInfo requests, that fails with the error 'Operation on mixed SRID geometries
A simple

select * from geometry_columns

returns that the SRID of pippo.geom is 0 (but should be 4326)

I tried also with this definition for the view

CREATE VIEW pippo AS SELECT st_setsrid(st_geometryn(shape,1), 4326) as geom FROM events

but I still have SRID=0 in the geometry_columns....any help?

Here is how the layer is defined in Geoserver:

Workaround

I found a workaround to fix not the problem itself, but to let Geoserver to perform spatial search on GetFeatureInfo requests:

instead of creating a view and publish it, I published an SQL-layer defined directly in geoserver. In the layer-definition form it allowed me to select the right gemoetry type and SRID and everything is working now!


texts: Interactively select parameters for snic segmentation in GEE

texts: Which functions/methods/ in GEE , one can use to play with or interactively choose snic segmentation parameters i.e., compactness,size,etc . I have tried ui.slider, linked map but need a code to follow please


texts: Creating Least Cost Path with cost distance analyse

texts: I'm trying to create Least cost Path analyse between two points in QGIS. I had DEM and from this I created Slope with various elevation values. Then I made two points with different layers. Both of them include only one point. When I want to run Least cost path algorithm, it appears problem starting point is no legal point, Then I changed the points but its same. Could someone help me with this issue?


texts:
A point layer with an ele attribute and another id_2 attribute. What I'm trying to do is to get the first 3 maximum elements of ele for each id_2

With a simple SQL query I can get the maximum element (just the max!) for each id_2:

SELECT geom, max(ele), id from layerA
group by id_2

How can I get the other 2 maximum elements?


texts:
I've been working on a bunch of layers of differing co-ordinates projections within a single dataframe.  Because I haven't really got my head around this, and because its all been working just fine, I didn't want to touch it.

However, one of the functions I was trying to perform wouldn't work (extract multi values to points) and I suspected this was because of coordinate conflicts.

My data frame properties show that I have layers with the following coordinates in it:

GCS_WGS_1984
NZGD_2000_New_Zealand_Transverse_Mercator
NZGD_2000_Transvers_Mercator
NZTM

I guess they should all be on the same.

I have just tried to change the projection of 3 layers (why didn't I check after the first one???) from GCS_WGS_1984 to NZGD_2000_New_Zealand_Transverse_mercator using "Define projection tool".  They've gone missing, even zooming to layer wont find it.

The properties of the missing layers show the extent in lat/long.  The properties of the layers I can see (haven't changed the projections of) show the extent in meters.

I'm wondering if this could be the issue, and if it is how to change it.  If its not this, what could it be?
Whats the best way forward?

I have re projected (as per advice) and can now see my layers on the same page - but they still don't overlap.  I suspect this is because of extents.  The attached images will show where the layers of extent type A and B respectively are showing.  Projecting has not changed this... is there a way to change it?  Maybe try another co-ordinate system?

]4


texts: I've got a TIFF imported as raster in my QGIS with KBS WGS84. Now I duplicate the layer and want to move that raster to the right edge of the first TIFF. So I just want to move / translate the TIFF to the right. Any suggestions how to do that?


texts: Buffering point shapefile using gdal/ogr: how to copy over input attributes to output shp?

texts:
Editing my previous question which asked the best way to buffer a point shapefile in python to add more detail and clarity to what I'm looking for.

I have found a method here to iterate through an input point feature class, buffer by a distance and write to an output polygon feature class:

def create_bufferShp_fromPointShp(inShpPath, bufferSize):

    buffShpPath = inShpPath.replace('.shp', '_buffer-{}m.shp'.format(bufferSize))

    inputds = ogr.Open(inShpPath)
    inputlyr = inputds.GetLayer()

    shpdriver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(buffShpPath):
        shpdriver.DeleteDataSource(buffShpPath)
    outputBufferds = shpdriver.CreateDataSource(buffShpPath)
    bufferlyr = outputBufferds.CreateLayer(buffShpPath, geom_type=ogr.wkbPolygon)
    featureDefn = bufferlyr.GetLayerDefn()

    for feature in inputlyr:
        ingeom = feature.GetGeometryRef()
        geomBuffer = ingeom.Buffer(bufferSize)

        outFeature = ogr.Feature(featureDefn)
        outFeature.SetGeometry(geomBuffer)
        bufferlyr.CreateFeature(outFeature)
        outFeature = None

    # Copy the input .prj file
    from shutil import copyfile
    copyfile(inShpPath.replace('.shp', '.prj'), buffShpPath.replace('.shp', '.prj'))

    return buffShpPath

The resulting polygon feature class is exactly what I'm looking for, but I would like to carry over the attributes from the input feature class to the output. These fields may change, so I am looking for an automated way. I added this portion to before the for loop:

    # Try to create new fields in the output shp
    for i in range(inputlyr.GetLayerDefn().GetFieldCount()):
        bufferlyr.CreateField(inputlyr.GetLayerDefn().GetFieldDefn(i))

This creates the fields, but I am unsure (outside of hardcoding the fields with outFeature.SetField()) of the best way to copy the input attributes of a feature to the output feature class.


texts: Performing Theil-Sen regression on raster data by R software

texts:
I have seven NDVI raster layers (in ERDAS img format) of an area. Each NDVI layer is the peak growing time of that year. I want to do trend analysis of NDVI values change over the years pixel by pixel. The seven time variables (years) are 1985, 1990, 1994, 1999, 2004, 2010, 2014. My intension is to use the Theil-Sen regression on the seven raster NDVI layers to get four raster outputs - slope, significance, adjusted significance and offset. Using the raster outputs, I will be able to generate a synthetic NDVI layer for any year in the time frame (for example 2007), which will help me for my further analysis.

Is there any way to get the four raster outputs (slope, sig, adj sig and offset) from the seven raster NDVI layers using Theil-Sen regression in R software?
If yes, what would be the potential package to be used?
What could be the potential code there for a less time consuming processing? I'm new in R.


texts: PostGIS dijkstra routing - specifying source and target as bigint instead of int4?

texts:
I'm currently working on a problem where I have to use PostGIS's pgr_dijkstra to calculate the shortest route from a starting point to a destination, and then use ST_AsGeoJSON to convert the route to a representation in GeoJSON. I came up with this query, which worked fine when I loaded a small number of geometries:

SELECT ST_AsGeoJSON(wkb_geometry) FROM (
SELECT seq, id1 AS node, id2 AS edge, route.cost, dt.wkb_geometry, dt.sidewalk_edge_id FROM pgr_dijkstra('
                SELECT sidewalk_edge_id AS id,
                         source::integer,
                         target::integer,
                         ST_Length(wkb_geometry)::double precision AS cost
                        FROM sidewalk_edge',
                469256344, 986535192, false, false) as route
                join sidewalk_edge  dt
                on route.id2 = dt.sidewalk_edge_id
) as routegeometries

sidewalk_edge is the table containing the geometries, wkb_geometry is the geometry column, and source and target are the columns containing source and column IDs (stored as varchar).

However, when I loaded the table with a larger number of geometries, some of which have bigger source and target IDs, I got this error:

ERROR:  value "5443376856" is out of range for type integer
********** Error **********

ERROR: value "5443376856" is out of range for type integer
SQL state: 22003

I tried changing source::integer and target::integer to source::bigint and target::bigint respectively, but PostGIS didn't like that:

ERROR:  Error, columns 'source', 'target' must be of type int4, 'cost' must be of type float8
********** Error **********

ERROR: Error, columns 'source', 'target' must be of type int4, 'cost' must be of type float8
SQL state: XX000

So it seems that my source and target IDs are to big to store as integers, but PostGIS will not accept anything other than integers. Is there any way to work around this problem?


texts: From routes of points make lines postgres ST_Makeline()

texts:
I have a dataset full of points that are routes. The data has time and id's with which I was able to ST_Union all the right points to form the multipoint for each route. Now I am trying to create a line of these multipoints. This is a harder challenge then I thought. ST_makeline just connects all points with each other and doesnt follow any order.

Now I think I need to reshape my multipoint to an array and provide this as an array to the ST_makeline(). But I am not sure if this is the right method.


texts: how to update values from python list into feature class

texts:
I have a feature class and python list. Feature class has 80 fields and 500 records. List also has 80 sub-list and 500 elements in each. So basically they are same dimension so to say.

I want to copy all records in python sub-list to corresponding fields in feature class. I was refering to an earlier post:

How to update attributes of a shapefile using list values?

I tried the code below but it doesnot work. there is some logic error can anyone correct me.

Code:

    fc = parameters[0].valueAsText   # feature class selected from python toolbox

    fname = parameters[1].valueAsText   #excel workbook selected from python toolbox.

    xl_workbook = xlrd.open_workbook(fname)
    sheet = xl_workbook.sheet_by_index(0)

    #here i am copying the excel sheet data to python list.
    l = [[sheet.cell_value(r,c) for r in range(sheet.nrows)] for c in range(sheet.ncols)]

    #sample of list :
    # l = [[u'OBJECTID', 1.0 , 2.0 , 3.0 , 4.0 ], [u'LENGTH', 56.29, 61.8 , 11.01 ,164.03]]

    #till here everything is working fine.
    #now i want to update the fields in feature class (fc) with the python list (l)

    i = 0
    cursor = arcpy.UpdateCursor(fc)
    for row in cursor:
        row[i] = l[i]
        cursor.updateRow(row)  #here i get 'valueError: Row: Invalid value for setting'
        i = i+1


texts: GeoExplorer provokes "401 Unauthorized" request and lacks login

texts:
Currently I have a login problem with my deployed GeoExplorer. Everything works fine in debug mode on my local VM with the geoserver from my remote server connected. Firebug gives the following output when loading GeoExplorer.

GET http://localhost:9080/geoserver/ows?SERVICE=WMS&REQUEST=GetCapabilities&TILED=true&VERSION=1.1.1 200 OK 140ms XMLHttpRequest.js (Zeile 220)

GET http://localhost:9080/geoserver/pdf/info.json 200 OK 18ms XMLHttpRequest.js (Zeile 220)

However, once I deploy GeoExplorer and move it to my remote server the login button on the upper right doesn't appear anymore. Instead the basic authorization window "GeoServer Realm" pops up and Firebug shows the following "401 Unauthorized" message:

GET http://Server_IP/geoserver/ows?SERVICE=WMS&REQUEST=GetCapabilities&TILED=true&VERSION=1.1.1 200 OK 227ms GeoExplorer.js (Zeile 396)

GET http://Server_IP/geoserver/pdf/info.json 200 OK 75ms GeoExplorer.js (Zeile 396)

GET http://Server_IP/geoserver/rest/imports 401 Unauthorized 202ms GeoExplorer.js (Zeile 396)

"NetworkError: 401 Unauthorized - http://Server_IP/geoserver/rest/imports"

The same problem applies to the GeoExplorer that was originally installed with the OpenGeoSuite, so it rather seems to stem from the server configuration instead of the deployment process.

I've tried various things in GeoServer security setup and found threads with similar symptoms, but couldn't find a solution for us yet. Maybe someone has another useful hint?!

Another idea: Can I configure the deployed GeoExplorer that it requests GeoServer over localhost instead of the server IP? Could this help?


texts: netCDF 4 file with relative time axis in arcpy

texts:
I'd like to process some netCDF 4 files with arcpy (ArcGIS 10.3.1). The original data comes from NCEP/NCAR project and contains wind vectors (uwnd/vwnd) for world-wide points (lon/lat) at certain time stamps.
The dimension "time" is given as relative time in

hours since 1800-01-01 00:00:0.0

If I have a look at the data with CDO (Climate Data Operators), the first four time stamps of the uwnd.10m.gauss.2016.nc are the following:

cdo sinfo uwnd.10m.gauss.2016.nc

2016-01-01 00:00:00
2016-01-01 06:00:00
2016-01-01 12:00:00
2016-01-01 18:00:00
...

Unfortunately arcpy has a problem understanding the relative time stamp and the following code

inNetCDF = r"\C:\temp\uwnd.10m.gauss.2016.nc"
ncFP = arcpy.NetCDFFileProperties(inNetCDF)
ncDim = ncFP.getDimensions()
for dim in ncDim:
    print "%s (%s)" % (dim, ncFP.getFieldType(dim))
    top = ncFP.getDimensionSize(dim)
    for i in range(0,top):
        print ncFP.getDimensionValue(dim,i)

delivers this as output for the dimension "time":

07.03.1800 17:50:24
07.03.1800 17:50:42
07.03.1800 17:51:00
07.03.1800 17:51:18
...

So instead of one time stamp every 6 hours ArcGIS calculates one time stamp every 18 seconds?!

It's not a big problem to recalculate the correct date within python, but I was wondering if there is a parameter I can set so that arcpy delivers directly the correct date and time. I've had the same problem with wrong date and time by using the MakeNetCDFFeatureLayer toolbox tool directly within ArcMap, and there I can't do any python calculations.


texts: Interpolation with Breaklines- How to Apply Breaklines as Surface Discontinuties

texts:
I can't seem to get Breaklines to behave as defined here with the QGIS 3 (3.10 today) built-in Interpolation (TIN or IDW) . In my case I am preparing a subsurface structure map with geologic faults (Breaklines) with vertical displacements. Contours should end at opposite sides of each fault reflecting the displacements.

I have seen one answer where the user moved TIN triangle vertices to end at Breaklines to control interpolation. This would be a tedious manual task for many subsurface maps. I am not a very proficient programmer so I am trying to use the built-in capabilities of QGIS. Most of the previous questions on Breaklines were submitted while the Interpolation Plugin was in use (prior to QGIS 3).

Is there some documentation how to prepare the Breakline Layer?

Are my expectations not consistent with QGIS developers' intentions?

I put together a simple map as a test. The following images are from that test. The GUI shows my selections for the interpolation. I tried the Use Z-coord selected as well.
Fault Layer Attribute Table

Resulting Map with contours.


texts: Is median center more influenced than mean center by value outlier in ArcMap?

texts:
I noticed something strange or seemingly paradoxical in ArcMap 10.2.2. The median center is less influenced by spatial outliers than mean center, this is a fact.
The thing is, I created 6 dispersed points and gave them similar attribute values. After this I calculated the mean and median center using those values as a weight field. The result was as I expected. The mean center point was influenced highly by the point with the highest value.

After this, I changed the value of a point and made it a high outlier. I recalculated the mean and median center and the result was very strange. The median seemed to be influenced more by that high outlier point. Why did this happen?

Thereafter, I tried the same thing but this time with an additional spatial outlier. At first I gave it a value similar to the others, and after this a very high value which is considered a high outlier value. The result was as I expected. The median center was less influenced by that outlier point (both with very high or similar to the others' value) than the mean center.

So the question is, does the spatial pattern of points play a role on weights? Does the nature of the outlier (spatial or value) point plays a role on the weights? Which of these exerts more influence to the position of the mean and median center points?

I'm uploading a screenshot which shows the first "experiment" I did.

The blue points are the means and the purple points are the medians. The heavily outlined points are unweighted. The smaller blue and purple points are the weighted with the high outlier value being under the selected point. The other two blue and purple points are weighted with the values of all the other points being similar

The reason for asking this question is because I do have empirical data uniformly spatially distributed with some values being outliers, and I'm trying to depict central tendency of several attributes. I noticed that "paradox" and didn't know which parameter to use. Mean center or median center?


texts: Reducing resolution of huge categorical raster in QGIS?

texts:
I'm not much of a GIS user, but I've got a big ascii file that I'm trying to convert into a

lat  long  data
x1   y1     d1
x1   y2     d2
x1   y3     d3
x2   y1     d4
...
etc

format, so that I can work with it in R.  I've got global national identifier data (i.e.: which country occupies what grid) at 2.5 arc-minute resolution, but I want to coarsen it into half-degree, because that is the resolution of the rest of my data.  I could do so in R if the file weren't so large -- its 3.9GB and I can't read it into memory.  A NetCDF file would obviously be ideal for me, but I can't find a netcdf version of the data I want.

The data is categorical -- country identifiers.  What I want to do is simple to take each 12x12 block of gridcells, evaluate that N=144 matrix, and return the country that appears the most times.  (I'd also like it to ignore NA values, unless all values are NA).  This would be easy enough in R if I could only open the darn text file.

So I opened QGIS, but there doesn't seem to be any obvious way to "coarsen" a raster.

Is there one?

Then it occured to me that I could convert the data to vector, and then convert it back to raster, specifying the resolution that I want.

BUT is there a way to do so that would allow me to specify that each gridcell should be .5x.5 degrees, going from -180 to 180 and -90 to 90?

This seems like a pretty gross hack.

Is there a better way?

I've decided that it makes more sense for me to just make the file into a netcdf using a remote machine with the memory for it.


texts: Aggregate census tract data to county level

texts: I have data on a census tract basis but I am having trouble aggregating them to county level. Is there any tool to do it in ArcGIS?


texts: How to generate a shapefile with UTM Grid precision level 100 m

texts: Does anybody know, if I can download from somewhere UTM Grid for zones 33 to 35 and U belt (for Poland) in shpfile or different format I can convert to shp. I need precision level 100 m. Is it possible to generate it in some GIS software?


texts:
I'm only just getting into dojo, so I have what most would probably consider a very basic question.  I've built my own legend in a div and added divs for the legend items.  This is all quite manual, but I just can't get my head around the dojo dijits.  I've written the code to toggle the divs on and off (just changing their background colour), but what I now need to add is the code to show or hide the layers from the map.

How can I use the ArcGISDynamicMapServiceLayer object I created to show and hide that entire service on the map after it's been constructed?  All the examples I've seen look at the service and then list the layers inside it.  I just want to add/remove or show/hide the entire service.


texts: Why pixel spacing measures as spatial resolution in Sentinel-1 SAR images

texts:
I have been processing Sentinel-1 SAR images-both medium and high resolution.
Looking into the documentation, the pixel spacing and resolution have been specified as:

High resolution level-1 GRD:

EW Swath -

resolution = 50x50 m
pixel spacing = 25x25 m

IW Swath -

resolution = 20x22 m
pixel spacing = 10x10 m

Medium resolution level-1 GRD:

EW Swath -

resolution = 93x87 m
pixel spacing = 40x40 m

IW Swath -

resolution = 88x87 m
pixel spacing = 40x40 m

However looking at the image, when I measure the length of the pixel, it is equal to the pixel spacing for both high and medium resolution images which is 10 m and 40 m respectively.

The images are shown below:

It is a level 1 GRD product. The high resolution image is in IW swath and the medium resolution in EW swath.

Is it then correct to say that the spatial resolution is 40 meters for medium resolution in EW swath, and 10 meters in IW swath?


texts: How do I repair NVIZ error: 'Can't find a usable init.tcl'?

texts:
How do I repair the following problem:

The error occurred while I was trying to run NVIZ in the GRASS environment in QGIS.  NVIZ however runs perfectly in GRASS.  On checking the directories listed, init.tcl is located there.


texts:
I am migrating some users from MapInfo to QGIS, and before they sign off on it, they want a very simple address search in the same way that we currently have one working in MapBasic.  I;ve been having a look around, but can;t seem to get my head around modifiying the quick find and find by attribute plugins which already exist.

Essentially they currently have the following:
A mapbasic routine which presents the user with portions of the address which they can fill in

Such as

field 1: Building name/number
field 2: Street Name
field 3: Locality
field 4: Town
field 5: Postcode

The process then queries a specific layer using like statements, returns a list of matching records which when clicked zoom to the relevant address.

They are adamant that using the QGIS query builder is going to be 'too complicated' for these users and want me to simplify the process in QGIS to a similar one to the above.

So we have a vector layer containing addresses available in the project - is it possible to mirror this MapInfo functionality - and if so, can anyone provide me with some guidance or point to a tutorial for doing this type of task.

Edit: 15:22 11/06/2013

To Clarify:
The addresses are held in a MapInfo tab file, which is loaded in to the project.  All of the fields that would be queried are of type string - there is an additional useful field called searchAddr which has field 2, 1 and 3 concatenated (split by commas).  The data also contains the coordinates as integers for each location.

So in terms of code (hopefully, if I use how it works in SQL to explain that will make it clear), we would be doing soemthing like:

Select * from <layer>
where searchAddr like '%<field1>%'
and searchAddr like'%<field2>%'
and searchAddr like '%<field3>%'
and Town like '%<field4>%'
and Postcode like '%<field5>%'

Therefore if the user puts in 42 as the House Number in field 1 and postcode AB12 3CD in field 5 the query would be

Select * from <layer>
where searchAddr like '%42%'
and searchAddr like'%%'
and searchAddr like '%%'
and Town like '%%'
and Postcode like '%AB12 3CD%'


texts: UpdateCursor is only updating the first letter from the value of my dictionary

texts:
I am using an update cursor to update an attribute field with the value of my dictionary.  Both the field and the value are strings.  I am confused why the update cursor is only updating the field with one character from the value of the dictionary. I am stuck and haven't been able to figure it out.
arcpy.env.workspace = r'\\lfkdshome\oit_home\Public\GIS\BWSO\GIS_UNIT\REQUESTS\Bioswales\RainGardentTesting\Bioswales_Archive.gdb'
mxd = arcpy.mapping.MapDocument("Current")
df = arcpy.mapping.ListDataFrames(mxd)[0]

BBL_Dict = {'247849': '3022000054', '664486': '4096010016', '924628': '4112140009'}

with arcpy.da.UpdateCursor("RG_164", ["BLDG_FID", "BLDG_BBL"]) as cur:
    for join in cur:
        if join[0] in BBL_DICT:
            join[1] = BBL_DICT[join[0]][0]
            cur.updateRow(join)


texts: Viz.json-based map will not show Line polygons

texts:
We are using a visualization and I importing it via a cartodb.createVisualization. Most of the tables work fine and the whole visualization works within the web console. However, we have one table which has a Line geometry. This geometry displays in the web console and appears in the website, but breaks when I use the .setSQL method. We have a large number of rows in this table (60k+). Often, the "loader" will hang. When I attempt to show this layer, the png GETs come back completely empty (transparent). Similarly, the grid.json requests show up empty, even when I have infowidnows I would expect to be populated with data. I am not getting any 404 errors.

When I request rows via the SQL API (like through sql.execute) the rows come back as expected. Also, I have tried using this visualization in my own website and in the example here, to make sure that it's not my mistake, but the problem is the same.

I don't know if it's related, but calling .show() on this sublayer also fails. Calling .show() on any of the sublayers will lead to unexpected behavior if I do it repeatedly (like instead of showing the correct sublayer, it shows a different sublayer that was hidden.)

Edit:  This geometry displays in the web console and appears in the website, but breaks when I use the .setSQL method. We have a large number of rows in this table (60k+ rows).

Update:  Cut the table from ~60k rows to ~14k and things are working as expected.


texts: Average Nearest Neighbor for X points instead of 1 point

texts: Is there a tool similar to the Average Nearest Neighbor tool where instead of calculating the average for the closest point it calculates the average for the closest X points?


texts: Im new to Linux, installed grass from the Ubuntu software centre. When I try to run it I get an error message saying there's an error with wxpython and to install grass gui. I have wxpython3, and setup the repositories from the grass download for Linux page. Then went to install grass-gui but it says multiple libraries are missing. Can anyone help?


texts: Opening images in TIFF format in QGIS 3.12

texts:
I am trying to open satellite imagery of several countries. The image descends in TIFF format. When I open it with QGIS 3.12 it opens in the form below. As far as I investigate, I do symbology settings, but it still doesn't happen. What should I do?

I am trying to open satellite imagery of several countries. The image descends in tiff format. When I open it with QGIS it opens in the form below.


texts: Dissolving specific polygons inside SearchCursor with Python in Arcmap

texts:
I have a layer with catchment polygons. For those that are smaller than a certain size, I want to merge with an adjoining polygon (namely the one downstream). In the attribute table I have a GridID, which i use to identify different polygons, and I can identify the downstream polygon by that attribute.

How can I, in this case inside a SearchCursor, merge these two polygons together?

I have tried with an UpdateCursor and arcpy.Dissolve_management() without getting it to work. Since Dissolve requires two tables with a value in common, I suppose I could just add a new field and make sure that this is the case for each iteration. But I can't use the same file as input and output (obviously), and I would like to get the file updated right away. If not, I might merge polygons that are of sufficient size since it has been merged with another earlier in the process, and after that meet the criterion.
In the case of an UpdateCursor, I think this is mainly for updating field values and not the geometry itself?

Is there a (more sufficient) method for this?

# inpgs = Catchment polygons shapefile

where1 = '"' + 'Shape_Area' + '" < 5000'
polyrows = arcpy.SearchCursor(inpgs, where1)

for prow in polyrows:
    grd1 = prow.GridID

    ...

    # grd1 = GridID on the current polygon
    # grd2 = GridID on the polygon downstream

    where4 = '"' + 'GridID' + '" = ' + str(grd2)
    polyrows2 = arcpy.UpdateCursor(inpgs, where4)

    for prow2 in polyrows2:
        arcpy.Dissolve_management(inpgs, outpgs, ???)
        inpgs = outpgs # ?

I use ArcInfo 10.0, Python 2.6 on Windows 7.

EDIT: I have now implemented the changes dmahr suggested, and made some more myself. That resulted in the following code:

# Create dictionary for all geometries
polyrows = arcpy.SearchCursor(inpgs)
geometryDictionary = {}

for prow in polyrows:
    geometryDictionary[prow.GridID] = prow.Shape
del prow, polyrows

# Create dictionary (nDD) which contains the GridIDs of the downstream catchment area
# as value and the current catchment area's GridID as key.

nextrows = arcpy.SearchCursor(indrln)
nextDownDictionary = {}
nDD = {}

for nrow in nextrows:
    nextd = nrow.NextDownID
    nextDownDictionary[nrow.GridId] = nrow.NextDownID

for k in dict.keys((nextDownDictionary)):
    clause = '"' + 'HydroID' + '" = ' + str(nextDownDictionary[k])
    nr = arcpy.SearchCursor(indrln, clause)
    for n in nr:
        nDD[k] = n.GridID
del n, nr

where1 = '"' + 'Shape_Area' + '" < 5000'
polyrows = arcpy.SearchCursor(inpgs, where1)
arcpy.Copy_management(inpgs, outpgs)

# Merge the small polygon into the downstream polygon
for prow in polyrows:
    grd1 = prow.GridID
    grd2 = nDD.get(grd1)

    if grd2:
        # grd1 = GridID on the current polygon
        # grd2 = GridID on the polygon downstream

        geometry1 = prow.Shape
        geometry2 = geometryDictionary[grd2]

        arcpy.Merge_management([geometry1, geometry2], tempMerged)
        arcpy.Dissolve_management(tempMerged, tempPgs)
        arcpy.Append_management(tempPgs, outpgs, "NO_TEST")
        grids.append(prow.GridID)

# After appending the dissolved polygons to the output file, the row for the merged polygon is being deleted
for g in grids:
    whereU = '"GridID" = ' + str(g)
    upd = arcpy.UpdateCursor(outpgs, whereU)
    for up in upd:
        upd.deleteRow(up)

del prow, polyrows
del up, upd


texts: Displaying specific areas specified by Global MapIt Area ID in OpenStreetMap using OpenLayers?

texts:
I've been customizing the Fix My Street platform for the two cities in my country, just to test everything out in it and how it works and maybe implement it for the whole country in the future. The two cities that will be covered are where we are going to do our testing. Right now, I'm already done with the minor customizations like the Admin boundaries, the areas where a person could only report and other front-end stuff. All of which was provided by the Fix My Street devs via their Cobrand module.

Seeing that they used OpenLayers to display OpenStreetMap, I've been searching in OpenLayers' mailing list and reading its documentation alot today, especially on how I could display the only two cities that my Fix My Street implementation covers. Kind of like what Fix My Street Brussels did with their map.

Hope you guys could give me a speed bump and point me to the right direction on where to start. I've already thought of other solutions but it uses other means to achieve this specific objective (e.g. Restrict user from dragging to areas not covered, display a grayscaled layer over my base layer except the areas that I covered.)


texts: Exceptions in arcpy using specific errors?

texts:
I have this code:

... for table in tableList:
...     with arcpy.da.SearchCursor(table, ["qname","tname","query"]) as cur:
...         for row in cur:
...             if row[2] not in ('',None):
...                 query = row[2]
...                 arcpy.Select_analysis(in_features=row[1],
...      out_feature_class=os.path.join(shapefile_output_folder,row[0]),where_clause=query)

If a certain file that meets at the loop doesn't exist I want to print the name of the file and continue without stopping.

The error when it doesn't find the file is :

  Runtime error  Traceback (most recent call last):   File "",
  line 8, in    File "c:\program files
  (x86)\arcgis\desktop10.6\arcpy\arcpy\analysis.py", line 90, in Select
  raise e ExecuteError: ERROR 000732: Input Features: Dataset
  TA2a does not exist or is not supported

So how should the exception be?
I tried:

try
 ... for table in tableList:
    ...     with arcpy.da.SearchCursor(table, ["qname","tname","query"]) as cur:
    ...         for row in cur:
    ...             if row[2] not in ('',None):
    ...                 query = row[2]
    ...                 arcpy.Select_analysis(in_features=row[1],
    ...      out_feature_class=os.path.join(shapefile_output_folder,row[0]),where_clause=query
)
except arcpy.RuntimeError:
print(table)

Can you suggest corrections?


texts: Interpolating climate data over DEM in ArcMap

texts:
I tried to interpolate climate data using IDW  by inserting DEM raster as a mask, but I got the Map of precipitation levels  without taking into account elevation data.
How can I get map of precipitation in according of DEM?
I am using ArcGIS Desktop 10.3.1.


texts: Offline map (Google/OpenLayers) for HTML5 (Angular) based hybrid application

texts:
I have a requirement where the customer has explicitly asked for offline support for the map. The rest of the application works well offline, but I am not sure how to get this done using Google or OpenLayers.

I have gone through here and I couldn't find any answer to what I was looking for. Please do note that I am only using angular.js and phonegap.js and nothing else.  I have worked with Google Maps as well as OpenLayers and I'm inclined to using these as I feel its much light weight and community support is awesome.

Q. Has anyone implemented or come across offline usage of Google Maps or OpenLayers. Other reliable map solutions are also welcome, but just need to be fast and mobile friendly.

What I have done:

I have downloaded the entire source repo from github (OpenLayers) and included in my application. Although it takes hell of a time loading in browser, the application when installed had no performance issue. The issue was even though I included all files in my app, OpenLayers still fetches the map images (tiles) from another server and it is like dynamic (I had asked related here: OpenLayers maps issue with SSL).

http://c.tile.openstreetmap.org/14/8745/5622.png
http://c.tile.openstreetmap.org/14/8745/5625.png
http://b.tile.openstreetmap.org/14/8745/5623.png

and so on.

I have not gone with Google Maps as everything I got is direct API call to Google server.


texts:
Let's say that I have two tables in PostGIS, tracks and cities, each table has two columns, in both cases the columns are called name and way, name is a string, way is a geometry.

I want to retrieve a list of tracks by name, something like:

SELECT name, way FROM tracks WHERE name = 'Mountain Trail';

But I also want to retrieve the name of the closest city to that track. Something like:

SELECT name, way, closest_city_name FROM tracks WHERE tracks.name = 'Mountain Trail';

But there is no closest_city_name column, I want to calculate the closest city as part of the query.

I'm familiar with ST_Distance and using joins and subqueries but I'm not sure how best to put it all together for this query.

Is there a way to JOIN cities to tracks using "closest"? Or should I query against the tracks table and a subquery on the cities table or vice versa?

Sorry this is quite a vague question, I'm just really looking for someone to point me in the right direction, perhaps an example query that you have used to do something similar?

Any suggestions greatly appreciated. Thanks!


texts: How to make input data matrix in flowmapper

texts:
I have problem with using flowmapper. I'm confused, how do I enter a text record? Can you give me an example of an appropriate text file?
I've been doing interpolation and made a contour of the ground water. How to create a flownet using flowmapper for my contour? This is my coordinate data:

Can you help me solve this problem?


texts: How to fix the grass7 package available from the ubuntugis-unstable PPA?

texts:
Out of the blue, grass stopped functioning this morning on one of the Ubuntu 14.04 systems I use:

$ grass
Starting GRASS GIS...
python: can't open file '/usr/lib/grass70/gui/wxpython/gis_set.py': [Errno 2] No such file or directory
Received EXIT message from GUI.
GRASS is not started. Bye.

After some fussing around I found that the grass meta-package available from the ubuntugis-unstable PPA was updated from 6.4.4 to 7.0.0 two days ago. Since I have this PPA registered in the sources.list file the grass meta-package got automatically upgraded. However, something went wrong, the grass-gui package is still version 6.4.4:

$ dpkg -l | grep grass
ii  grass                                                 7.0.0-1~exp2~trusty                                 all          Geographic Resources Analysis Support System (GRASS GIS)
ii  grass-core                                            7.0.0-1~exp2~trusty                                 amd64        GRASS GIS core components
ii  grass-doc                                             7.0.0-1~exp2~trusty                                 all          GRASS GIS user documentation
ii  grass-gui                                             6.4.4-1~trusty5                                     amd64        GRASS GIS graphical user interfaces
ii  libqgisgrass2.8.0                                     1:2.8.0+20trusty                                    amd64        QGIS - shared grass library
ii  libqgisgrass2.8.1                                     1:2.8.1+20trusty                                    amd64        QGIS - shared grass library
ii  qgis-plugin-grass                                     1:2.8.1+20trusty                                    amd64        GRASS plugin for QGIS
ii  qgis-plugin-grass-common                              1:2.8.1+20trusty                                    all          GRASS plugin for QGIS - architecture-independent data

As it happens, the new grass-gui is marking a number of dependencies that can not be installed:

$ sudo apt-get build-dep grass-gui
Reading package lists... Done
Building dependency tree
Reading state information... Done
Picking 'grass' as source package instead of 'grass-gui'
Note, selecting 'libtiff5-dev' instead of 'libtiff-dev'
The following packages have unmet dependencies.
 libcairo2-dev : Depends: libcairo2 (= 1.13.0~20140204-0ubuntu1) but 1.13.0~20140204-0ubuntu1.1 is to be installed
                 Depends: libcairo-gobject2 (= 1.13.0~20140204-0ubuntu1) but 1.13.0~20140204-0ubuntu1.1 is to be installed
                 Depends: libfontconfig1-dev (>= 2.2.95) but it is not going to be installed
                 Depends: libglib2.0-dev but it is not going to be installed
 libproj-dev : Depends: libproj0 (= 4.8.0-2ubuntu2) but 4.8.0-4~saucy2 is to be installed
E: Build-dependencies for grass-gui could not be satisfied.

Is there any hope of getting these dependencies installed? Do I need some extra PPA for it? Or should I simply revert back to version 6.4.4?


texts: Calculating distance between points following flow direction on ArcMap?

texts:
I delineated flow directions in a watershed on ArcMap (10m DEM). Then I want to calculate distance between two points 'following' the flow direction in each raster cell, not a straight line distance.

Is there any way to do it automatically on ArcMap?

Added description:
The below image shows an example distance (black line) that I want to calculate. Those arrows represent a flow direction in each cell (I changed the symbology to show those directions by arrows). The distance is not on a channel.


texts: QGIS Python Script - Routing along Waypoints

texts: I used Anita Graser's script point_layer_to_route.py which worked like a charm until I updated to 2.6. With QGIS 2.6 I always get objects without any features. Can anybody tell me if he/she had the same issues?


texts: Using update cursor with list of shapefiles?

texts:
I'm still pretty new to using cursors. What I'd like to do is create a list of shapefiles. Within those shapefiles I would like to populate a field "NAME" with the shapefile name less the .shp extension. I am using 10.1 and my code is below:

bufferCalc = arcpy.ListFeatureClasses("*.shp")
for calc in bufferCalc:
NOTE: This space contains calculations for other fields
with arcpy.da.UpdateCursor(calc, ("PROJNAME")) as cursor:
    for shp in arcpy.ListFeatureClasses():
        cursor.updateRow(os.path.splitext(shp))

The other calculations I need run smoothly however, the error I get when the cursor is ran states "StopIteration: iteration not started". This has got to be something simple I'm doing wrong but can't seem to figure it out.


texts:
I'm using QGIS to render some air quality modelling results.  I'm having issues with the OpenLayers Plugin whenever I have another vector layer active.  The plugin works if it's the only layer in the project, but displays the following error whenever I add another layer.  Has anyone else experienced this or know what to do?

An error has occured while executing Python code:
Traceback (most recent call last):
  File "C:/Users//.qgis//python/plugins\openlayers_plugin\openlayers_ovwidget.py", line 250, in __signal_webViewMap_loadFinished
    self.__refreshMapOL()
  File "C:/Users//.qgis//python/plugins\openlayers_plugin\openlayers_ovwidget.py", line 265, in __refreshMapOL
    action = "map.setCenter(new OpenLayers.LonLat(%f, %f));" % self.__getCenterLongLat2OL()
  File "C:/Users//.qgis//python/plugins\openlayers_plugin\openlayers_ovwidget.py", line 275, in __getCenterLongLat2OL
    pntCenter = coodTrans.transform(pntCenter, core.QgsCoordinateTransform.ForwardTransform)
QgsCsException: forward transform of
(11049.424326, 84358.680211)
PROJ.4: +proj=longlat +datum=WGS84 +no_defs +to +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs
Error: latitude or longitude exceeded limits

Python version:
2.7.2 (default, Jun 12 2011, 15:08:59) [MSC v.1500 32 bit (Intel)]

QGIS version:
1.8.0-Lisboa Lisboa, 6416f38

Python path: ['C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python', 'C:/Users//.qgis//python', 'C:/Users/hamishh/.qgis//python/plugins', 'C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python/plugins', 'C:\\Windows\\system32\\python27.zip', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\DLLs', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\plat-win', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\lib-tk', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\qgis\\bin', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\site-packages', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\site-packages\\PIL', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\site-packages\\win32', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\site-packages\\win32\\lib', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\site-packages\\Pythonwin', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\Python27\\lib\\site-packages\\wx-2.8-msw-unicode', 'C:\\Users\\\\.qgis\\python\\plugins\\QuickMultiAttributeEdit/forms', '~/.qgis/python', '.', 'C:\\PROGRA~2\\Quantum GIS Lisboa\\apps\\qgis\\python\\plugins\\fTools\\tools']

Sorry for the long snippet, thanks for your help.


texts:
I'm trying to get a list of the recent earthquakes, with latitude and longitude from this API: https://apps.pdc.org/msf/rest/services/global/pdc_active_hazards/MapServer/6

I can't rely on a server-side library as I have to use only JavaScript in an HTML page.

My problem is that I can't find the list of the events browsing the API


texts: Calling marker-size from GeoJSON file in leaflet.js

texts:
I found marker-size: medium in the properties section in the GeoJSON file. so I want to call that marker-size, to be used in the PHP Leaflet code.
So I need a solution for using this marker size and how can I solve it easily
===========geojson data=========================

     {
            "type":"FeatureCollection",
            "features":[
                {
                    "type":"Feature",
                    "properties":{"markercolor":"#1e7e7e","marker-size":"medium","marker-symbol":"","title":"Indonesia"},
                    "geometry":{"type":"Point","coordinates":[107.9296875,-6.664607562172573]}
                },
                {
                    "type":"Feature",
                    "properties":{},
                    "geometry":{"type":"Point","coordinates":[114.873046875,4.565473550710278]}
                },

========leaflet php code===========

pointToLayer: function(feature, latlng) {
   var smallIcon = L.icon({
                      iconSize:feature.properties.marker_size,
                      iconAnchor: [13, 27],
                      popupAnchor:  [1, -24],

   });

   return L.marker(latlng, {icon: smallIcon});
},
        onEachFeature:onEachFeature
    }).addTo(mymap);
    });


texts: How to Update the EPSG Database of Gdal on Ubuntu?

texts:
I have a Ubuntu 16.04 Server, with GDAL 2.1.3, released 2017/20/01.

I need to project a Raster to EPSG:7755 which was added way back in 2016.

When I use this EPSG code with gdalwarp, I get the following error:

ERROR 6: EPSG PCS/GCS code 7755 not found in EPSG support files.  Is this a valid
EPSG coordinate system?
ERROR 1: Translating source or target SRS failed:
EPSG:7755

How do I Update the EPSG Database of Gdal?


texts:
I'm trying to add a vector layer from vector source , but getting error like this

AppComponent.html:1 ERROR TypeError: this.adaptAttributions_ is not a function
    at Object._ol_source_Source_ (source.js:41)

here is my code ,

import {Component, OnInit} from '@angular/core';

import OlMap from 'ol/map';
import OlXYZ from 'ol/source/xyz';
import OlTileLayer from 'ol/layer/tile';
import OlView from 'ol/view';
import OlProj from 'ol/proj';
import OlProjection from 'ol/proj/projection'
import OlExtent from 'ol/extent';
import OlTileGrid from 'ol/tilegrid/tilegrid'
import OlLayer from 'ol/layer/layer';
import OlLayerImage from 'ol/layer/image';
import OlSource from 'ol/source/source';
import OlSourceOSM from 'ol/source/osm';
import OlTileWMS from 'ol/source/tilewms';
import OlVector from 'ol/source/vector';
import OlLayerVector from 'ol/layer/vector';
import OlFormatGeoJSON from 'ol/format/geojson';
import OlLoadingStrategy from 'ol/loadingstrategy';
import OlImageWMS from 'ol/source/imagewms';
import OlStyle from 'ol/style'
import Proj4 from 'proj4';
//import * as ol from 'ol';
@Component({
  selector: 'app-map',
  templateUrl: './map.component.html',
  styleUrls: ['./map.component.css']
})

export class MapComponent implements OnInit {

  map: OlMap;
  source: OlXYZ;
  layer: OlTileLayer;
  view: OlView;

  constructor() {
  }

  ngOnInit() {

      var vectorSource = OlVector({
        format: new OlFormatGeoJSON(),
        url: function(extent) {
    return 'http://localhost:8080/geoserver/workspace/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=workspace:roadcentreline&outputFormat=application/json&srsname=EPSG:32640&'+
    'bbox=' + extent.join(',') + ',EPSG:32640';

        },
        strategy: OlLoadingStrategy.bbox
      });

      Proj4.defs("EPSG:32640","+proj=utm +zone=40 +datum=WGS84 +units=m +no_defs");
      var proj27700 = OlProjection.get('EPSG:32640');

       var vector = new OlLayerVector({
        source: vectorSource,
        style: new OlStyle.Style({
          stroke: new OlStyle.Stroke({
            color: 'rgba(0, 0, 255, 1.0)',
            width: 2
          })
        })
      });

      this.view= new OlView({
                    center: [376581.6867, 2679756.7678],
                    projection: 'EPSG:32640',
                    maxZoom: 19,
                    zoom: 15
        });

        this.map= new OlMap({
        layers: [vector],
        target: document.getElementById('map'),
        view: this.view
      });

  }
}

i used same version of openlayers(4.6.5) used here

working here, but not for me.


texts: Reducing coordinate precision in MSSQL WKT using C#?

texts:
When you use the AsText() function on geometry objects in MSSQL Server 2008+ you get something like this:

"LINESTRING (8.8446177236171337 47.721989356913951, 8.844617919281065 47.721977134398223, 8.844589303417802 47.721976942728766, 8.84454742773029 47.7219801783058, 8.8445167302356325 47.721989634438941, ...)"

The precision of the coordinates is far to great to be useful (for me). For wgs84 coordinates the seventh decimal place is worth ~11 mm. After that, it gets more or less meaningless.

For my workflow I need to put a lot of WKT's in a Textfile. To reduce file size I would like to reduce the coordinate precision to six decimal places.

I could process the textfile and reduce precision in my backend (c#, entity framework) but that seems like a rather costly task.

Is there a way in MSSQL Server to create WKT from geometry objects with a defined number of decimal places?


texts: I'm looking for a good technique to take points, which mark the locations of fence post holes, from a Google map layer to then precisely mark with paint where to drill the holes. I noticed that my neighbor's fence isn't straight, and I want mine to be geographically straight. :-) I'm looking for a simple device or iPhone app to find the points I have for the fence boundaries out in the field.


texts:
I am currently working on the effect of land use on streams in very small catchments (<30km²) in Germany. We have accurate data (ATKIS) on general land use patterns (forests, agriculture, urban areas etc.). We also try to incorporate information about crop types (CAPRI data) in the agricultural areas into our analysis. However the CAPRI data we use is quite coarse: Aggregated between 2001 and 2003 and a spatial solution of 1000x1000m. http://esdac.jrc.ec.europa.eu/content/european-food-safety-authority-efsa-data-persam-software-tool

Is there any open source data on specific crop types in Europe (or Germany) for certain years other than CAPRI?


texts: Using OFFSETB in viewshed analysis of ArcGIS Desktop?

texts:
I am confused by the description of OFFSETB in wiewshed analysis (ArcGIS 10.1).

The help documentation says:

  The OFFSETB item indicates a vertical distance in surface units to be
  added to the z-value of each cell as it is considered for visibility.

Still, to me is rather obscure.

I understood that OFFSETA should be the value representing the height of the observer point: that is, if an observer is 2 meters high, and (s)he is on a spot 10 meters above sea level, the total elevation of the observer point will be 12m.

As to OFFSETB: is that the height of a target object?
What this actually represents?
Or is it the height (in the landscape) below which the viewshed is not calculated?


texts: Symbolizing layer with graduated colors by a field in Excel File

texts:
I have a field in Excel file called percentage, I want the governorates colored by the percentage value, so the more value is the more red become darker in the governorate.

So how can I do it in ArcMap?


texts: How to get street view using OpenLayers with Bing map

texts:
I am looking into solution for displaying street view using OpenLayers with Bing map. I couldn't find any examples on OpenLayers website. They have examples which displays satellite map using different imaginary set.
Is it possible to use OpenLayers to display street view using Bing Map? If yes then how.


texts: problems starting publishing services with ArcGIS for Server 10.2

texts:
I have just installed ArcGIS for Server 10.2 on my server and the publishing services will not start. Has anyone had issues with this before? The error message that I receive after the service fails to start is:

  Service failed to start  System.PublishingTools.GPServer: java.rmi.UnmarshalException: Error unmarshaling return header; nested
  exception is: java.net.SocketException: Connection reset

If anyone has some insights it would be much appreciated.


texts: Count max number of pixels for image in GEE and then calculte numbers from this number

texts:
I have a code that creates image collection and then calculate what is the biggest number of pixels in image (e.g the image with maximum spatial data).
I managed to calculte this number and to use it to create new image collection with relevant images.
I want to create range of numebr of pixels that are needed in an image in order to join this image collection. I decided that I can "accept" images that are at least 75% size of the biggest image.
I have a code for that, but whenever I try to compute max_number.multiply(0.75) I get error:

  max_number.multiply is not a function

Those are therelevant parts from my code:

var dataset = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
                  .filterDate('2013-01-01', '2019-12-02')
                  .select('B4', 'B5','pixel_qa')
                  .filterBounds(geometry)
                  .map(maskL8sr);

//clip the dataset according to the geometry
var clippedCol=dataset.map(function(im){
  return im.clip(geometry);
});

// Get the number of images.
var count = clippedCol.size();
print('Count: ',count);

//here I do calculations -clip image collection, add NDVI so my new
//my new dataset called now "withNDVI"

//Filter according to number of pixels

var ndviWithCount = withNDVI.map(function(image){
  var countpixels = ee.Number(image.reduceRegion({
  reducer: ee.Reducer.count(),
  geometry: geometry,
  crs: 'EPSG:4326',
  scale: 30,
  }).get('NDVI'));

  return image.set('count', countpixels);
});

print('NDVIwithCount:',ndviWithCount);

var max = ndviWithCount.reduceColumns(ee.Reducer.max(),  ["count"]);
print(max.get('max'));
var max_number=(max.get('max'));
var almost_max=max_number.multiply(0.75);

My end goal is to be able to compute numbers with my var max_number so in the end my script knows automatically for every polygon to choose only the correct images


texts:
I have 14 layers, which are polygons. They represent the channel of a river at different periods. I need to fill the gaps that exist once I overlay the polygons, but I can't change the properties of the original layers because I need to reclassify the resultant polygon on the base of the colour in order to do more analysis. I have tried to merge all 14 polygons, but the new polygon doesn't have the same properties of each individual polygon.

Do you know how I could fill the gaps between the 14 polygons?

FINAL SOLUTION:

Finally I have done it by creating two layers; one by union all layers without gaps, and other by union the layers with gaps. Then, I used Symmetrical differences tool and I obtained the layer with gap polygons.


texts: Erase within a feature class - Planar enforcement

texts:
I´m using ArcGIS v10 ArcInfo.

I´m working with a feature class that has 5 features which are waterdepths (1 to 5m). These features are overlapping each other. What I need, is that they are not overlapping. To achieve this I used two differend methods that are leading to the same goal:

I extracted each feature in a single feature class (1, 2, 3, 4, 5). 1 is the biggest overlapping 2, 3, 4 and 5. So I erased 1 with 2, 2 with 3, 4 with 5, merged them and had the feature class I needed.
I ran the union-tool (geoprocessing) on the feature class itself. The problem is, that all overlapping parts are created additionally as new features. Instead of the 5 features I need, I get 10 features more which I have to sort and erase by hand.

My question: Is there a tool or method that does this without so many steps by hand? I´m able working with QGIS and ArcGIS v10 ArcINFO.


texts: Why does ST_Area in SQL give a different value to $area in QGIS?

texts:
The area value I'm getting using the SQL term 'ST_Area' is markedly different from that gained from QGIS functions like $area, the area measure tool and area as a derived feature (all of which agree with one another).

I've tried three different CRS, one geographic and two projected, project CRS matches the sqlite layer. Changing the features shows that areas are definitely being updated.

Could anyone explain what 'ST_Area' is measuring and why it's not the same as Qgis measurements?


texts: Selected Feature not highlighted in map of CMV (Configurable Map viewer)

texts: I have copied the whole code of Find widget and added it into another widget. First widget is working fine and highlighting the feature(creating feature from geometry), second widget is also working fine. i.e. it is returning me the accurate and proper number of results, and same zoom extent but it is not drawing the feature, which are return by the query.


texts: I want to print the results of the identify features, I know that there is a way to do it with the Python console, but the thing is that i have to set up an easy approach for one who is not into QGIS or Python. Has anyone an idea or an easy  workaround?


texts: Converting .map() from JavaScript API to Python API?

texts:
I wrote this code using the Earth Engine JavaScript API. The code clips the image collection using variable "table," which is a feature collection of about 12,000 points. I get the length of the image collection and use that number to create a a variable "sequence." I use this variable to create a nested map loop of sorts. I grab a single image and get a list of coordinates from that image. Then for each set of coordinates in the list, I store a bunch of information from that point into a dictionary and create a feature out of it. I need to convert this code to Python, but I am not sure how to convert the lines of code that use the map() function twice ("var getImgCoords and var getMetadata"). The variable "sequence2" changes for each loop, so I thought I could use the iterate() function instead of map().

var TerraClimate = ee.ImageCollection("IDAHO_EPSCOR/TERRACLIMATE"),
var table = ee.FeatureCollection("users/adrianom/CAN_AK_Coordinates")

var SearchArea = table;
var SearchAreaGeom = SearchArea.geometry();
var terraclimate_bounds = TerraClimate.filterBounds(SearchAreaGeom);

  // Clip collection to search area, more precisely than filterBounds above.
  var icClip = terraclimate_bounds.map(function(img){
    var i = img.clip(SearchAreaGeom);
    return(i);
  });

  var icClip_len = icClip.size();
  var col_list = icClip.toList(icClip_len);

  var sequence = ee.List.sequence(0, icClip_len.subtract(1));

  var getImgCoords = sequence.map(function(i){
    var img = ee.Image(col_list.get(i));
    var coords = img.geometry().coordinates();
    var coords_len = coords.size();
    var sequence2 = ee.List.sequence(0, coords_len.subtract(1));

    var getMetadata = sequence2.map(function(j){
        var coor = coords.get(j);
        var geom = ee.Geometry.Point(coor);
        var reduce = img.reduceRegion(ee.Reducer.mean(), geom, 4638.312116386398);

        var dictionary = ee.Dictionary({
          x: ee.List(coor).get(0),
          y: ee.List(coor).get(1),
          Start_Date: ee.Date(img.get('system:time_start')),
          End_Date: ee.Date(img.get('system:time_end'))
        });

        var combine_dict = reduce.combine(dictionary);
        var feature = ee.Feature(geom, combine_dict);
        return(feature);
    });

    return(getMetadata);

  });


texts:
I am trying to plot a shapefile on top of a raster using the command readOGR. My current code is:
    CEZ <- readOGR(".", "Shapefile/CEZ1")

This is followed by the error code

Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv = use_iconv,        :
Cannot open data source

I am unsure what this code means - Is it a problem with the shapefile, the command or the way I have structured the path?


texts: Using arcpy from the command line on ArcGIS Linux Server 10.1

texts:
We installed ArcGIS SERVER (not Desktop!) 10.1 on a RHEL (Red Hat) machine.  ArcGIS Sever is supported for Linux.

However, it appears they don't make it easy to script in python, from the command line, making use of the arcpy module.

There are many things awkward about this.  The arcpy library that they use is actually for Windows, with some external patching. In fact, it comes with its own implementation of python. So pathnames look a little weird. You have to be the arcgis user to even get access to python or arcpy.   So here is how my code looks:

import arcpy
"""
From http://resources.arcgis.com/en/help/main/10.1/index.html#//0155000005mm000000

The Python libraries installed by an ArcGIS Linux Server installation is Windows 64-bit Python. Therefor, stand-alone Python scripts should always use the Win\
dows path separator (\) when accessing data and other Python modules. Also, be aware that there is a Z:\ which is mapped to the Linux systems' /. Within a Pyt\
hon script you can access data and other Python modules using Z:\ mapping. For example, data in /myuser/myproject/data can be accessed as z:\myuser\myproject\\
data.
This does not apply to Python scripts you publish from ArcGIS for Desktop because the publishing process converts all paths to the correct format.

"""
"""

I'm using:

sudo su arcgis

/opt/arcgis/server/tools/python ./arcpytest.py

"""

dataDir=r"z:home\arcpy-tmp"+"\\"
print dataDir

shp='tgr39105lkA.shp'
arcpy.FeatureToLine_management(dataDir+shp, 'temp_links')

What I get back is:

$ /opt/arcgis/server/tools/python ./arcpytest-1.py
z:home\arcpy-tmp\
Traceback (most recent call last):
  File "./arcpytest-1.py", line 28, in <module>
    arcpy.FeatureToLine_management(dataDir+shp, 'temp_links')
  File "C:\Program Files\ArcGIS\Server\arcpy\arcpy\management.py", line 2343, in FeatureToLine
    raise e
arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Input Features: Dataset z:home\arcpy-tmp\tgr39105lkA.shp does not exist or is not supported
Failed to execute (FeatureToLine).

Has anyone tried to use arcpy in this direct way? It's a REALLY powerful server, so we put ArcGIS server on it, but want to use ArcPy for programming.


texts: Choosing projection to georeference historical map to?

texts:
I'm using ArcGIS.

I have a historical map from 1884. It has bathymetry information which I would like to compare to modern bathymetry models with as much accuracy as possible. I don't know what projection to use to retain the best accuracy for comparison.

I've tried georeferencing it once using a WGS_84 projection, and again using a NAD_83 projection for that zone (since this is a small area).Naturally, the map warps differently depending on which projection I've used.

A problem with this map as well, is that I can't space the GCPs very evenly across the raster because the entire bottom half of the map is water. So while I can get the top half georeferenced to fit nicely (using the roads, lighthouses and even some natural features), I'm afraid the bathymetry information will be incorrectly distorted. That's why I really want to choose the best projection.

Which one do you think I should use?

Do you have any ideas on how I can keep this most accurate?


texts: How to project the result of a QueryTask in ArcGIS Javascript API?

texts: I have created a web map that displays point features based on a given date range. Our current projection is Oregon State Plane South. My boss wants me to change the map to use a Web Mercator projection. I am loading our web mercator basemap fine, but I'm having trouble projecting and displaying the points. What is the most efficient way of doing this?


texts:
I have some points I would like to map out in 3D. The points were gotten from a GPS module athat was attached to a Raspberry Pi.
This is the settings I used before starting qgis2threejs:

This was the output it gave me:

The purple vector building's settings are as shown:

The issue is that the altitude values for the pink spheres were taken on ground level. So why is it floating in the air like that?
Does the CRS matter? I enabled on the fly CRS transformation already. The land that I am mapping out is relatively flat, around 20-40m elevation from sea level.
I am not understanding this problem! Desperate for help!
UPDATE
I found the problem!
The problem is that the Z coordinate mode adds up the height
"relative to DEM" + "Altitude" + Addend
Is it possible to remove the "relative to DEM" to the equation?


texts: Stable QGIS version running on Windows 7

texts:
I am having trouble with QGIS 3.10 constantly crashing on a Windows 7 machine.
We are looking at previous versions to install to get it to work. Upgrading to Windows 10 isn't an option.
What is the latest version that works on Windows 7?


texts: Is it possible to retrieve thematic layer information for a Layer in Mapbasic?

texts:
I found this code

 windowinfo(FrontWindow(), WIN_INFO_CLONEWINDOW)

This function gives information (mapbasic) about the front Map window which can be used to regenerate the map.

  Is there any function in mapbasic 10.0 function that gives code that can be used to generate a
  thematic layer.


texts:
I am very new to arcgis in general, and am having trouble stiching several rasters together. I know that the mosaic tool can do this, but when I try to create a new mosaic dataset I get the error message "you do not have the necessary license to execute the selected tool".

The spatial anaylist extension is turned ON, and I have checked it several times, and it is otherwise working. Is there something else besides the spatial analyist I need?

I couldn't figure out how to post on Esri's forum, and hoped someone here could help me. I am getting desperate! Any help is appreciated. Thank you!


texts: Copy only Minimum value which is minimum within group in QGIS?

texts:
In my data attribute:

'FID_Net' is 'unique value/subset' and 'class' is the 'value'. And in my new data sheet, the values of column 'class' is not repeating.
Now I have to copy that unique value 'FID_Net' to new field 'New_FID' of only one row which value 'class' is minimum within group 'FID_Net', others row should be 0 or NULL.

As shown in below picture, how can i do in Qgis 3.x?


texts:
I'm trying to create a point IFeatureClass, but my code doesn't work.

Public Function CreatePointFeatureClass(ByVal featureClassName As String, _
                                         ByVal featureWorkspace As IFeatureWorkspace) As IFeatureClass

    ' Instantiate a feature class description to get the required fields.
    Dim fcDescription As IFeatureClassDescription = New FeatureClassDescriptionClass()
    Dim ocDescription As IObjectClassDescription = CType(fcDescription, IObjectClassDescription)

    Dim fields As IFields = New FieldsClass
    Dim fieldsEdit As IFieldsEdit = CType(fields, IFieldsEdit)
    Dim fieldUserDefined As IField = New FieldClass
    Dim fieldEdit As IFieldEdit = CType(fieldUserDefined, IFieldEdit)

    ' Set up geometry definition for the Shape field.
    ' You do not have to set the spatial reference, as it is inherited from the feature dataset.
    Dim geometryDef As IGeometryDef = New GeometryDefClass
    Dim geometryDefEdit As IGeometryDefEdit = CType(geometryDef, IGeometryDefEdit)

    ' By setting the grid size to 0, you are allowing ArcGIS to determine the appropriate grid sizes for the feature class.
    ' If in a personal geodatabase, the grid size is 1,000. If in a file or ArcSDE geodatabase, the grid size
    ' is based on the initial loading or inserting of features.
    geometryDefEdit.GeometryType_2 = ESRI.ArcGIS.Geometry.esriGeometryType.esriGeometryPoint

    ' Set standard field properties.
    fieldEdit.Name_2 = "SHAPE"
    fieldEdit.Type_2 = esriFieldType.esriFieldTypeGeometry
    fieldEdit.GeometryDef_2 = geometryDefEdit
    fieldEdit.IsNullable_2 = True
    fieldEdit.Required_2 = True
    fieldsEdit.AddField(fieldUserDefined)

    ' Add a X Double field to the required fields.
    Dim field1 As IField = New FieldClass()
    Dim fieldEdit1 As IFieldEdit = CType(field1, IFieldEdit)
    fieldEdit1.Name_2 = "X"
    fieldEdit1.Type_2 = esriFieldType.esriFieldTypeDouble
    fieldsEdit.AddField(field1)

    ' Add a Y Double field to the required fields.
    Dim field2 As IField = New FieldClass()
    Dim fieldEdit2 As IFieldEdit = CType(field2, IFieldEdit)
    fieldEdit2.Name_2 = "Y"
    fieldEdit2.Type_2 = esriFieldType.esriFieldTypeDouble
    fieldsEdit.AddField(field2)

    ' Add a Elevation Double field to the required fields.
    Dim field3 As IField = New FieldClass()
    Dim fieldEdit3 As IFieldEdit = CType(field3, IFieldEdit)
    fieldEdit3.Name_2 = "Elevation"
    fieldEdit3.Type_2 = esriFieldType.esriFieldTypeDouble
    fieldsEdit.AddField(field3)

    ' Use IFieldChecker to create a validated fields collection.
    Dim fieldChecker As IFieldChecker = New FieldCheckerClass()
    Dim enumFieldError As IEnumFieldError = Nothing
    Dim validatedFields As IFields = Nothing
    fieldChecker.ValidateWorkspace = CType(featureWorkspace, IWorkspace)
    fieldChecker.Validate(fields, enumFieldError, validatedFields)

    ' The enumFieldError enumerator can be inspected at this point to determine
    ' which fields were modified during validation.
    ' Create the feature class.
    Dim featureClass As IFeatureClass = featureWorkspace.CreateFeatureClass(featureClassName, validatedFields, _
                                    ocDescription.InstanceCLSID, ocDescription.ClassExtensionCLSID, esriFeatureType.esriFTSimple, _
                                    "SHAPE", "")

    Return featureClass

End Function


texts:
I've got a Layer where old data versions are still in the layer. We've got the field "name" and the field "Version". The problem is that not every Name gets a new Version at the same Time. So some names are still version 1, others are version 5.
Now I have to find the newest Version for every Name and put only them in the new Layer.
Right now it takes me two hours to do so manually (in Excel CSV)...
Each feature has one entry for "name" and "Version".
So names are repeating: several features have the same name, but different versions.
And now I want to compare all features with the same name to see which one has the newest version.
I am not deep into QGIS and using a German version but I'll try to explain: There is a Vector function which extracts by attribute; that is not working since it only asks for one attribute I thought about using the field calculator but I can't work out how to do such complex things like asking other IDs if they've got the same name. I have no other ideas yet.


texts: My WMS Server returna a template without image instead of map's image

texts:
I have a server on localhost (using apache and mapserver in windows). I want to create a wrapper and do it using this document : http://mapserver.org/cgi/wrapper.html - Apache ReWrite rules (using Apache mod_rewrite)

I create a .htaccess file and write this code:

RewriteEngine on
RewriteRule map?(.*) cgi-bin/mapserv.exe?map=C:\\OSGeo4W\\apache\\htdocs\\rasht\\RashtMap.map&$1

My WMS Server load with this address correctly and return an image of map:

http://localhost/cgi-bin/mapserv.exe?map=C:\OSGeo4W\apache\htdocs\rasht\RashtMap.map&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&LAYERS=Base,Street&STYLES=&CRS=EPSG:32639&BBOX=373204.15,4128154.94,377063.25,4131409.17&WIDTH=800&HEIGHT=500&FORMAT=image/png

But with wrapper address only load my template without any image.I want a response like previous address. I type this address as wrapper address :

http://localhost/map?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&LAYERS=Base,Street&STYLES=&CRS=EPSG:32639&BBOX=373204.15,4128154.94,377063.25,4131409.17&WIDTH=800&HEIGHT=500&FORMAT=image/png

What is wrong?


texts: Join 2 layers based on overlapping features in QGIS

texts:
I have two layers: (1) zipcodes (2) census tracts

I am analyzing a selection of 100 census tracts. Is there a way to process which census tracts are in each zipcodes based on their overlap?

i.e. for zip code 21223 there is census tracts 250203, 250204 250301 etc.

There are many census tracts (red border) in each zipcode (black border).

But each census tract (purple with red border) could also be in many zipcodes (yellow layer).

Is there a way to generate this zipcode-census tract list or join/create attribute tables based on the two layer's overlap? The zipcode layer has a column of 'zipcode' and the census tract has a column of 'census tract id,' but there is no mutual information between the two attribute tables.


texts: Error for using a buffer to insert features only when the Feature Class's domain is set

texts:
I have constructed the Spatial Reference and Coordinate system for a new FeatureClass.  I have set the extent of the FeatureClass and use the same spatial reference for features when I add them to the FeatureClass.

I set the FeatureClass domain by:

ISpatialReference pSpatRef;
ISpatialReferenceFactory pSpatialFactory = new SpatialReferenceEnvironment();
IGeographicCoordinateSystem pCoord =
            pSpatialFactory.CreateGeographicCoordinateSystem((int)esriSRGeoCSType.esriSRGeoCS_WGS1984);
pCoord.SetDomain(Math.Floor(XMin), Math.Ceiling(XMax), Math.Floor(YMin),
                Math.Ceiling(YMax));
pSpatRef = pCoord as ISpatialReference;

I have verified that the Domain is set properly by looking at the Feature class's properties.

I would like to insert the new features using IFeatureCursor with IFeatureBuffer.

My code for doing this is:

 IFeatureCursor pCursor = pGPXFeatureClass.Insert(true);
 IFeatureBuffer pBuffer = pGPXFeatureClass.CreateFeatureBuffer();

 foreach (string fileName in fileList)
 {
     GPXtoPolyline parseGPX = new GPXtoPolyline();
     parseGPX.GPXParser(fileName, pSpatRef);

     ///find the index of the name and shape attributes
     int intNameIndex = pGPXFeatureClass.FindField("Name");
     int intShapeIndex = pGPXFeatureClass.FindField("Shape");

     pBuffer.set_Value(intNameIndex, parseGPX.strName);
     pCursor.InsertFeature(pBuffer);
 }
 pCursor.Flush();

Here is my problem.  When I run my procedure without using the SetDomain method it works without fault.  However, when I use the SetDomain method, I get the following error:

  An unhandled exception of type 'System.Runtime.InteropServices.COMException' occurred in GEOG489_FinalProject_PGW.dll

  Additional information: TabletPC inking error code. Queue is full (Exception from HRESULT: 0x80040238)

This error is triggered by the command:

pBuffer.Shape = parseGPX.pPolylineGeometry;

What is the cause of this error?


texts: PyQGIS gives TypeError: QgsExpression.prepare(): argument 1 has unexpected type 'QgsFields'?

texts:
This question is very related to the unanswered Retrieve features from QgsFieldExpressionWidget pyqgis

I have a layer and I've build a QgsExpression that involves a field of that layer. What I'm trying to do is to evaluate the expression and get the results as a single list.

The problem is that I'm trying to do this with QGIS master (so the future 3 version) and the instruction of the PyQGIS cookbook are outdated because some API are changed.

Just to have an example:

vl = iface.activeLayer()
exp = QgsExpression('PH > 7') # where PH is the name of a field
exp.prepare(vl.pendingFields())

# then I get this error

Traceback (most recent call last):
  File "/usr/lib/python3.5/code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
TypeError: QgsExpression.prepare(): argument 1 has unexpected type 'QgsFields'

Someone has already successfully used the QgsExpression for QGIS master?


texts: Does joining CSV file to shapefile overcomes the 255 character limitation of shapefile?

texts:
I have a CSV file in which some fields contain values greater than 255 characters. I want to ask that wether joining the CSV file with a shapefile overcomes the 255 character limitation of shapefile and displays the data in CSV completely.
If not, then what is the solution for it?


texts: How to always see all band values with Identify Features of QGIS?

texts:
Does anyone know if it is possible to get the Identify Features tool in QGIS to show the band values for all the layers shown in the results window?

To clarify, whilst I have several bands "showing" in the results window, the actual values for all but the topmost layer are folded away.

Having to unfold the info each time I click on a new location is getting pretty tiresome.


texts: PostGIS: Faster way to ST_Clip() rasters with a polygon?

texts:
I want to clip an image with a polygon from postGIS rasters.
What I do:
 - Import raster into postGIS(raster2pgsql, tile-size: 500x500, each tile contains 12 bands landsat 8)
 - Import a shapefile into postGIS with SRID = 4326.
I have 21309 tiles in my raster table, and 60 pylygons in my shape table. Now, I want to clip a polygon:

SELECT oid,
   lowrite(lo_open(oid, 131072), png) AS num_bytes
FROM (
      VALUES (lo_create(0),
              ST_AsPNG(
                         (SELECT ST_Clip(ST_Transform(
                                           (
                                           SELECT ST_Union(foo.rast)
                                            FROM
                                              (SELECT rast
                                               FROM ldcm15, vnm_adm2
                                               WHERE gid=18
                                                 AND ST_Intersects(rast_geom_4326, geom)) AS foo),4326), vnm_adm2.geom,true)
                          FROM vnm_adm2
                          WHERE gid = 18)))) AS v(oid,png);

rast_geom_4326 is where I stored rast's polygon. Because, I don't want to use ST_Transform(rast, 4326) in the query, it's pain to wait for re-projection.
However, it still takes a long time to get the result.

I'm using PostGIS 2.1 on PostgresSQL 9.3.


texts: Selecting polygons that are to left (in west) of point in ArcMap?

texts:
I am working with ArcMap 10.3 and I have to select polygons that are left (in the west) of a point layer. I want to try it with "select by location" but I am not sure if that is really possible.

Is there a way that I can handle this in ArcMap?


texts: Find if lat/long point is inside an area

texts:
I have got data using OpenStreetMaps Ways of many areas.  I would now like to write a function which would return "true or false" if an single latitude/longitude point is found in that area or not.

I have already seen many answers to this question such as this answer however none of these deal with the issue when the Greenwich meridian line is crossed.  Most of my area searches will occur in London which means I need a function that can cope with both negative and positive values for the lat/long points in the area.

Is there a known way on how to solve this?


texts: Search and concatenate Landsat bands with Python

texts:
I have a number of Landsat geotiff (*.TIF) file in a folder. The problem is simple, I need to search certain Landsat bands from that list and make a composite tif image. I am using Orfeo Toolbox's CLI utility. So far I have this

import os
dir = "C:\Documents and Settings\Landsat_image"
def fileList():
    for root, dirnames, filenames in os.walk(dir):
        for tifFile in filenames:
            for i in range(3):
                if tifFile.find('B' + str(i)) != -1 and tifFile.endswith('.TIF'):
                    inputImages = os.path.join(dir, tifFile)
                    outputImages = os.path.join(root, tifFile) + '_concat.tif'
                    process = os.system(
                    'otbcli_ConcatenateImages -il %s -out %s' % (inputImages, outputImages))
print fileList()

The otbcli_ConcatenateImages needs the inputs in following format

otbcli_ConcatenateImages -il GomaAvant.png GomaApres.png -out otbConcatenateImages.tif

But maybe it is getting multiple entries for out.

ERROR: Invalid number of value for: "out", invalid number of values 4

I know the problem is silly but I am still a noob.


texts:
I've been using QGis for quite a few years. I normally update QGis to the latest version once it is released on qgis.org, without uninstalling previous versions. I've never encountered any problems doing this.

I recently installed 2.8.1 and it was working fine. I thought I'd fiddle with a bit of customization of QGis by way of switching off certain toolbars, menus, setting a default project to open, etc,  for someone to use QGis as a simplified viewing tool for gis data. I decided to try this on one of the older versions (chugiak). It seemed to work ok while in the program, and even while restarting the program. At some point, not sure whe, possibly on a reboot,  or  when I tried using 2.8.1 again I got an error (screenshot attached).

Since then I've tried everything I can think of. Reinstalling, uninstalling all versions and reinstalling 2.8.1. Still didn't work. I even tried deleting the windows registry variables pertaining to QGis (I may not have done this completely correctly) and again reinstalling, yet I still get the same error.

I get the splash screen, and I can see this screen saying things like "setting up gui", "loading plugins", etc. I don't get to the tips screen and the program gui doesn't come up.


texts:
I have a point table I have sorted by mmsi (it is the identifier) and by datetime columns.

  mmsi numeric,
  x1 double precision,
  y1 double precision,
  datetime timestamp without time zone,
  geom geometry(Point,3035)
  linetime2 (interval)

For each mmsi I would like to create a set of segments every 2 consecutive rows/points, skipping segments which have more then 6 hours of difference between the starting and the ending point AND more then 30 kilometers of distance between the starting and the ending points. the data are as follow:

Concerning the maximum distance between two consectutive points of the same mmsi I could calculate the line length in a second step and then delete lines longer than 30 km.

Concerning the maximun time interval I have calculted the time difference (column linetime2) between two consecutive points/rows using a lag function:

select mmsi,
       datetime,
       datetime - lag(datetime) over (order by mmsi, datetime) as linetime
from test_st_12_v1;

with cte as (select mmsi, x1, y1, datetime, datetime - lag(datetime) over (partition by mmsi order by mmsi, datetime) as linetime
from test_st_12_v1)
update test_st_12_v1 lt
set linetime2 = ct.linetime
from cte ct
where lt.mmsi = ct.mmsi
and lt.datetime= ct.datetime;

Now, how can I create lines every two consecutive points with the same id (mmsi) skipping segments with a linetime2 > 6 hours? E.g. in the example I dont' want a segnemt between the points 9 and 10, 15 and 16, 83 and 84, and of course between 84 and 85 because I have another id number (mmsi).
Any idea?


texts: Importing Camera images to ArcGIS or QGIS with correct location, scale and oritentation

texts:
I have several very large image datasets from a AUV survey, essentially aerial imagery, that I need to bring into a GIS environment and have them display at the proper spatial scale, orientation and location.  The images had no GPS EXIF metadata embedded, but with the help of EXIFTool I have appended this information.  X,Y,height above seafloor(Altitude), vehicle heading, focal length.

We have a MATLAB script that will mosaic raw images together, but this survey was poorly conducted and there are gaps in the along track overlap and no cross-track overlap.

My hope is that I can bring them into ArcGIS as a mosaic_dataset or something similar in QGIS.  They do not need to be stitched together, just displayed in their proper location, orientation and spatial scale.

How do I define the spatial scale of an image (Hopefully batch for the ~8000/dataset), and the orientation? Each image should cover real world area of 2.5m x 4m.


texts: ArcGIS Field Calculater DateAdd() gives Expected end of statement error

texts:
I need to add seconds in an numeric column(TS_sec) to a datetime column(TS) in ArcGIS field calculator. I am using Vbscript.
But i am getting following error:
ERROR 999999: Error executing function. Expected end of statement Failed to execute (Calculate Field (16))

I am using following VBscript code:

Code block
----------------------------------------------------
dim TS_seconds as date
if [TS_sec] = Null then
TS_seconds=Null
else
TS_seconds=DateAdd ("s", [TS_sec] , [TS] )
end if
----------------------------------------------------
Expression ETA_TS=
TS_seconds


texts: Activate / deactivate direction arrows in Mapillary JS

texts:
I need to activate / deactivate direction arrows on Mapillary images in a web application built on Mapillary JS.

I've tried to build a little code sample to explain ....

Here you are the code

    <meta charset='utf-8' />
    <title></title>
    <meta name='viewport' content='initial-scale=1,maximum-scale=1,user-scalable=no' />

    <link href='https://unpkg.com/mapillary-js@2.5.0/dist/mapillary.min.css' rel='stylesheet' />

        Select option:
        <select onchange="onSelect.call(this, event)">
          <option value="with">With arrows</option>
          <option value="without">Without arrows</option>
        </select>

The initial situation is WITHOUT the direction arrows as you can see because I start with this configuration

           {
                component: {
                    cover: false,
                    direction: false
                }
            }

I'd like to modify it using a selection above the image and if I select "With  arrows" I'd like to show the arrows ....

I've tried to manage this situation with this code ...

              case "With arrows":
                 var theComponent = mly.getComponent('direction');
                 theComponent.configure({direction: true});
                 break;

that execute without errors but no changes are on my Mapillary image ....

Suggestions / examples?


texts: I'm a newbie using QGIS 2.12 and I'm trying to split a line in QGIS and place points every x meters. There is a tool in ArcGIS (divide by length add in). How to achieve this in QGIS?


texts:
I'm using 10.0 because that's what our customer that I'm attempting to build a custom add-in for is using.

I am following the "Walkthrough: Building custom UI elements using add-ins" from the ArcObjects help file in Visual Studio 2010 C#.

I'm attempting to install the build on a virtual that has ArcGIS engine runtime 10.0 using the esriAddIn.

I get an error "Installation failed.  Matching ArcGIS installation not found.".

The walkthrough references ESRI.ArcGIS.Desktop.AddIns.  I obviously don't have desktop on my virtual so in VS2010, I set the dll reference for that to "Copy Local" and copied it over with the esriAddIn installation package.

Is it possible to use the Desktop.AddIns reference and install it on a machine that only has runtime and not desktop?


texts: Is there a way of changing reference scale of existing dimension feature class?

texts:
I have an existing dimension feature class that I would like to change a reference scale of. Is that possible at all?
Or do I have to always create a new one, with desired reference scale, and then import styles from existing one?


texts: Creating metadata style in ArcGIS Desktop from GeoNetwork plugin schema?

texts:
ArcGIS 10 does not have a  style metadata that meets my specification.

I have a GeoNetwork plugin schema for metadata generation as the Geospatial Metadata Profile of Brazil.

Is it possible to generate a metadata style in ArcGIS from it?


texts: Combine duplicate rows but keep class information?

texts:
I'm stuck trying to figure out how to avoid combining and deleting duplicate values manually. I'm working with parcel data, and based on some property information, each parcel is either in group a, b, c, multiple groups, or no groups at all. I had to do some select by location processing to get the lists of the parcels for each group.

I want to create a master list that combines the records from each group, removes duplicates, but still retains the group information. For example, if I had:

John Doe, 123 Main Street, Group A
John Doe, 123 Main Street, Group B
John Doe, 123 Main Street, Group C
Mary Smith, 321 Pleasant Street, Group B
Mary Smith, 321 Pleasant Street, Group C
Alex Jones, 111 Park Ave, Group C

I would want the output to be:

John Doe, 123 Main Street, Group A, Group B, Group C
Mary Smith, 321 Pleasant Street, Group B, Group C
Alex Jones, 111 Park Ave, Group C

I already created GROUP_1, GROUP_2, and GROUP_3 field in the respective attribute tables and have tried merge, union, append, and even some manual select by location add to selection and cannot figure it out for the life of me.

I am using ArcMap.


texts: Geoserver unable to execute a PostGIS stored procedure

texts:
I've a PostGIS stored procedure that does an insertion in a table, then it returns a value. Here it's:

BEGIN
    INSERT INTO cities (id, name) VALUES (17, 'koko');
    RETURN 2;
END;

I was able to call the stored procedure from geoserver using "Configure new SQL view..." the view retuns the value correctly (which is 2), but no insertion is done in the DB.

What is the problem?


texts: How to use expression in QgsVectorlayer label by using python?

texts:
how to show the two field values as a label on features. I have done it for a field called "ErrorDescr" with this i need to show another field values. Other field is "MarkerId"

layer.setCustomProperty("labeling", "pal")
    layer.setCustomProperty("labeling/enabled", "true")
    layer.setCustomProperty("labeling/fontFamily", "Arial")
    layer.setCustomProperty("labeling/fontSize", "10")
    layer.setCustomProperty("labeling/fieldName", "ErrorDescr")
    layer.setCustomProperty("labeling/placement", "2")

is it possible to bring the output like this:


texts: Is there a free dataset of historical geocoordinates of ancient cities and places?

texts:
Is there a freely available dataset (csv, geojson, topojson, etc.) with geocoordinates of ancient cities, for example of the Roman Empire? Ideally with additional information like population size?

I've searched GeoCommons and NaturalEarth but found nothing.


texts: Getting wrong values of Flow direction in ArcGIS?

texts:
I am trying to generate flow direction from DEM. However I am getting some vales as 80 and 129, though it should be only from 1, 2, 4, 8,16, 32, 64, 128

What could be the reason for this error?
How to resolve it?

Edit: Original DEM is at 90m resolution, whereas I am generating fdr at 0.25 degree spatial resolution by setting the cell size as 0.25 degree in raster analysis.


texts: STATISTICS not allowed in the Attribute Table

texts:
I created a shapefile, calculated geometry then went to grab the stats, the statistics was ghosted, not allowing me to grab those numbers.  I am using 10.1 ArcView.

I deleted the field and created a new field in ArcCatalog but no luck, same results.  Not sure why that is happening.

Some help, please.


texts: How to create a buffer zone for notification

texts:
I am working on a project, for this I need to find out abutter notification that generates a list of properties within the buffer zone from the selected property with user input distance in radius and generate a mailing list of each property displayed by the buffer result.

We are using GeoServer server and PostGIS, and need to do this on the server side coding and display result in OpenLayers viewer. Can anyone guide me with PostGIS spatial query to find  parcels from the selected parcels?

Distance measured for this is from the parcel outline not from the center of property.

I am new to PostGIS.


texts: Using Calculate Value to Replace UTM83 to NAD83 in file name within Model Builder

texts:
I am looking to use the Calculate Value tool to replace the UTM83.sid input aerial file name with an output aerial file name using the new projected coordinate system of NAD83.jp2

What do I use for an expression to get the desired result?
Is the replace expression applicable?


texts: I'm a new user to QGIS, trying to solve a problem that I've encountered in ArcGIS Online. I want to run a 'Find Nearest' analysis, however the number of points I'm routing from is greater than 1,000. Is there a simple way to undertaken this task in QGIS? I need the distances to be based on the road network as opposed to crow fly distance.