texts: Selecting max(cursor) Value returns incorrect max from arcpy.da.UpdateCursor

texts:
I'm working through some da.UpdateCursor scripts to make updates to a fgdb. Although the script below works successfully as a rule, it returns the incorrect max value from one feature class.

The feature class in question has been rebuilt numerous times but I still get the second highest value returned not the highest. Note this occurs in this feature class on the three columns with values (OBJECTID, Shape_Length, Shape_Area), the only other column is Shape.

Can anyone provide some guidance please as I have researched this and can't get an answer after some hours of experimenting across numerous other feature classes.

My only thought at this point, 'is there a ceiling on max?' given that  min(cursor) returns the lowest value correctly on the three columns.

Using Win 7, ArcGIS 10.2, Python 2.7.5 (yes, I'm in a mixed OS environment).

    import arcpy, sys, os
    from os import sep

    ws = env.workspace

    with arcpy.da.UpdateCursor(ws +sep+ 'Library_362', "Shape_Area") as cursor:
        for row in cursor:
            highValue = max(cursor)
            print highValue

    # returns [7693891.489877664]
    # s/be   [10054353.128549]


texts: How to create randomly points within polygons for each row of a dataframe matching a Polygon ID with R

texts:
My question is similar to this post, but because there are no "perfect" answer and above all no reproducible example I open this question.

I have a "SpatialPolygonsDataFrame" and a data frame. Both have a variable of common factors.
I would like to generate a random point for each observation of the data frame within the according polygon.

I am starting to learn GIS with R. I read the first chapters of Bivand et al. 2013 but i am still strugling with the basics. I tried with the function spsample from the library(sp) but I am puzzling with the variable n of this function: it's an (approximate) sample size

In my example there are 3 polygons and 30 observations. The common variable is MatchID. I need a SpatialPolygonsDataFrame with 15 points in the first polygons ("abcd"), 10 in the second ("efgh") in 5 in the third ("ijkl"). See

table(df$MatchID)

Example:

 library(sp)
# Definition of the CRS
poly.crs <- CRS("+proj=utm +zone=36 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")

# Definition of 3 Polygons
poly.a <- Polygon((matrix(c(-4.653724,1.2210259,-3.803160, 1.3799638, -3.245480, 0.1665987,-4.666098, 0.1097523, -4.653724, 1.2210259), nrow=5, ncol=2, byrow=T)))
poly.b <- Polygon((matrix(c(-5.820343, 2.675320,-5.427519, 3.062975,-4.701119, 2.819967,-4.555540, 1.855489, -5.050758, 1.388308, -5.783673, 1.572882,-5.820343, 2.675320), nrow=7, ncol=2, byrow=T)))
poly.c <- Polygon((matrix(c(-3.758639, 2.873654, -3.273958, 3.417311, -2.213099, 2.935320, -1.972031, 1.945992, -3.033510, 1.279312, -3.709359, 1.844633, -3.758639, 2.873654), nrow=7, ncol=2,     byrow=T)))

# Making a SpatialPolygons
polys.a = Polygons(list(poly.a), "pa")
polys.b = Polygons(list(poly.b), "pb")
polys.c = Polygons(list(poly.c), "pc")
Spolys = SpatialPolygons(list(polys.a,polys.b,polys.c), 1:3,  proj4string=poly.crs)

# Making a SpatialPolygonsDataFrame
data.Spolys<- (data.frame(MatchID=c("abcd","efgh","ijkl"), row.names=row.names(Spolys)))
Poly <- SpatialPolygonsDataFrame(Spolys, data.Spolys, match.ID = TRUE)

## Data frame that should be converted in a SpatialPointsDataFrame thanks to the variable "MatchID"
df <- data.frame(
          MatchID=c("abcd","abcd","abcd","efgh","efgh","ijkl"),
          V1 = 1:30,
          V2 = "a"
)

I need a function to generalize the idea below (looking for the matching with the factor) and giving exactly the number of observation I need (and not an approximate random size)

spdf.a <- SpatialPointsDataFrame(spsample(Poly[Poly$MatchID=="abcd",], n = nrow(df[df$MatchID=="abcd",]), "stratified"), df[df$MatchID=="abcd",])


texts: How to add multiple text elements in ArcMap using ArcPy

texts:
I've been able to add a dynamic text element (with an element name of 'TestName') that reflects a layer name and that changes for each corresponding exported layer to a GIF. I am now trying to add a second text element (with an element name of 'TestName2') for a map title but am running into issues. Below is the code I am working with. CMD prompt produces this error:

C:\UserS\JameS>export5.py
[<map layer u’Div1_Irrig_1987_2’>, <map layer u’Div1_Irrig_1987_2’>, <map layer
u’Div1_Irrig_1956_0’>, <map layer u’citylyr’>, <map layer u’boundlyr’>]
Traceback (most recent call last):
   File "C:\Python27\ArcGIS10.2\Scripts\TestScripts\export5.py”, line 42, in <module>
   elm.text : ""
   File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\arcobjects\7base.py", line 89, in _set
      return setattr(self._arc_object, attr_name, cval(val))
RuntimeError: TextElementObject: Error in setting text

C:\Users\James>

for lyr in lyrlist:
        lyr.visible = False
        #if lyr.name in layers and lyr.name <> layer:
        #   lyr.visible = False
        if lyr.name == "citylyr":
            lyr.visible = True
        if lyr.name == "boundlyr":
            lyr.visible = True
        if lyr.name == name:
            lyr.visible = True
            layerOn = lyr.name

            for elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
                if elm.name == 'TestName':
                    elm.text = lyr.name
                if elm.name == 'TestName2':
                    elm.text = 'Example Map Title Name'
                else:
                    elm.text = ""

    if layerOn:
        arcpy.mapping.ExportToGIF(mxd,r"{}\{}.gif".format(GIFPath,layerOn))


texts: Using the model builder in ArcGIS, how do I iterate through rasters in a model?

texts:
Right now I'm attempting to iterate through rasters collected in a file geodatabase, in a model that I've created using the model builder in ArcGIS. I don't have the ModelBuilder Toolbox, otherwise this would be a lot simpler to do.

Is there a way (without using the Iterate Rasters tool) that I can loop through the rasters in a gdb?

Thanks,
 Renee


texts: How to get parks in American Samoa from GeoNames

texts:
I have this query which does not return anything:

http://api.geonames.org/search?featureCode=PRK&username=genadinik&country=US&style=full&adminCode1=AS&maxRows=1000

Should I change it?  Is American Samoa considered its own country?  Un-american Samoa? lol.

Thanks!


texts: Batch processing DTM with Hillshade in QGIS 3

texts:
I'm using the Hillshade tool from GDAL library [Raster analysis] in QGIS 3.4.

I can't use the batch process with creative options like compress=LZW or profile=Geotif. It doesn't work.

The tool works only an image for time. Has someone found the same problem?


texts: Layer overlay problem in Openlayers /Geoserver

texts:
Hallo I am trying to create a Web Map servise using Openlayers and Geoserver. I want Open Street Map (OSM) as back ground and that both vagar_alla and adresspt layers to be tranparent, so that it will be possible both layers at the same time. I have the following code sofar. The problem is I OSM layer appears only at the begining zoom-in it disappears AND the top most layer obscures the other layers beneath. Any help?

<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/html4/loose.dtd">

        Lund Map

location

                --> Click on the map to get feature info


texts: Keeping the original file name when batch processing GDAL translate?

texts:
I would like to convert .tif files to Idrisi or ENVI format. I would like to do a batch process, keeping the original filename. When I choose: "run as a batch process", as I can see, I still have to set the output name and directory for each file?

Is there a faster way to do this?


texts: How to create a radius search in Leaflet

texts:
My objective is to create a radius search tool using the Leaflet API. It will behave like How many POIs within the radius from a given point?

Any tutorial or any blog forum link will be helpful to me.


texts: How to get custom colormap onto 3D NVIZ in QGIS GRASS

texts:
I am trying to create a 3D version of my 2D raster using QGIS GRASS NVIZ. I have been following these steps: http://linfiniti.com/2010/12/3d-visualisation-and-dem-creation-in-qgis-with-the-grass-plugin/

I have no problems creating the 3D model, but my colors always get lost.  This is frustrating because the second option on the NVIZ module asks for a raster to use for color.  I've tried multiple times to get NVIZ to create the 3D model with the color of the raster indicated and I always get the Pseudocolor or grayscale colormap.  The link provided in the above directions (http://linfiniti.com/2010/12/a-workflow-for-creating-beautiful-relief-shaded-dems-using-gdal/) claims to explain how to map colors, but does so only through command line procedures.  I am not yet familiar enough with GRASS command line to take that route and meet my deadline.

This GISSE question did not really explain how to do it through the NVIZ module:
How to change the raster color range in GRASS NVIZ?

Any suggestions would be much appreciated.

-mike


texts: Django raster, not add raster to BD Postgis

texts:
I am new in this world, and I am trying to make Django Raster work, I have already followed all the instructions of:
https://django-raster.readthedocs.io/en/stable/installation.html
But at the moment of creating the raster with the administrator, (I attached a .tiff image) I receive the following error::

  The process does not have access to the file because it is being used
  by another process.

Also when attaching a tiff file a new tiff file with the same name is generated plus a randomized code example: test.tiff generates a test_1dcv.tiff

I am currently working with GeoDjango + Postgis in Window and Celery is False in settings.


texts: Setting transparency in an SLD file for a 4 channel GEOTIFF

texts:
I'm struggling to set the a background colour to transparent in an SLD file for a 4-channel GEOTIFF using GeoServer 2.1.4.

I think I've set my background colour to transparent alpha layer in a graphics package correctly.

I've set up an SLD file to read Red/Green/Blue channels of my GEOTIFF = this resolves colours correctly for the RGB image = but still shows the background black colour (see below).

When I do a GetFeatureInfo request on the WMS, I can see an ALPHA_BAND=0 where I want the transparency and ALPHA_BAND =255 where I do not want transparency (see attached figure).

I was hopeing to do something like setting the ALPHA_BAND=0 to Opacity=0 in the SLD or something. However, there's no  option in SLD schema that I can see ?

I've read the SLD cookbook and the SLD reference documentation.  Alpha transparency is referred to but only with ColorMap' property. Can I mix and match 'ColorMap' and ChannelSelection' properties in the SLD ?

Am I going about this the right way ?  Thanks in advance for assistance!

This is the SLD I'm using:

    <?xml version="1.0" encoding="ISO-8859-1"?>
    <StyledLayerDescriptor version="1.0.0" xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc"
      xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
      xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.0.0/StyledLayerDescriptor.xsd">
      <NamedLayer>
            <Name>bathy example</Name>
            <UserStyle>
              <Name>RGB image 3 channels</Name>
              <Title>RGB image</Title>
              <FeatureTypeStyle>
                    <Rule>
                      <RasterSymbolizer>
                            <ChannelSelection>
                              <RedChannel>
                              <SourceChannelName>1</SourceChannelName>
                              </RedChannel>
                              <GreenChannel>
                              <SourceChannelName>2</SourceChannelName>
                              </GreenChannel>
                              <BlueChannel>
                              <SourceChannelName>3</SourceChannelName>
                              </BlueChannel>
                            </ChannelSelection>
                      </RasterSymbolizer>
                    </Rule>
              </FeatureTypeStyle>
            </UserStyle>
      </NamedLayer>
    </StyledLayerDescriptor>


texts: Calculating slope in percentage with QGIS

texts:
I'm working with QGIS, v. 2.0.1. I have several soil types. They are (amongst others) described with slope indications. I already calculated the slope from the DTM in degrees, but as description of the soil types is in percentage, I need to calculate the slope in percentage.
I tried two things:

I tried the analysis tool from the raster section, but it didn't work. The result was just black with values from 0 to 0.
The same happened when I tried to calculate the percentage in the raster calculator but the result was the same.

Maybe I need to say that the terrain is pretty plain.
What else can I do?
Update
I tried the raster calculator and did what @radouxju suggested. It gave me a black picture with values from 1756.64 to 1756.64. But as soon as I click in the picture the values are "normal", say 5.88 or 13.224 or something like that. The same happens when I try the gdaldem tool but with values from 0 to 0. I also tried to change the contrast enhancement (like stretch to min/max) but it didn't change anything.
Update 2
I just found out that maybe I have to set the scale to 111120  but that doesn't work neither (result: black, values 0-0).


texts: Insert a "scale box" of given side length in the QGIS layout?

texts:
QGIS allows displaying a scale bar, indicating distances. I would like to complement it with a box showing a given area, e.g. 1 hectare. How can I display this in a layout, ideally as an extension of the scale bar?
Example:
.----------------.

|........1ha.......|

'-----------------'    100   -----------------    200   -----------------    300 -----------------      400    -----------------   500 m


texts: Error in commit changes to Spatialite database

texts:
I have developed a cadastral application with QGIS 2.18.13. The cadastral plots are stored as polygons (WKB type "Polygons") in a Spatialite database, a number of attributes is included.

This worked well for a while but now I am facing a commit error when I try to save changes?

Nothing is stored in the database. The error message is as follows (translated from Dutch):

Commit of changes has failed for layer Erf

Errors: ERROR: 1 item not added.

 Provider errors:
  SQLite error: unknown cause
SQL: INSERT INTO "Erf"("geometrie","pkuid","Perceelnummer","Naam","Voornamen","Beroepen",
"Woonplaatsen","Artikel kadastrale legger","Soort van de eigendommen",
"Nummer der gebouwde eigendommen","Inhoud perceel","Bijeengetrokken soort van bebouwing",
"Voorlopige klassering","Polder- of Dijklasten",
"Tarief ongeb Klasse 1","Tarief ongeb Klasse 2",
"Tarief ongeb Klasse 3","Tarief ongeb Klasse 4","Tarief ongeb Klasse 5",
"Kad ink geb eigendommen","Kad ink bijeengetrokken","Aftrek wegens polder- of dijklasten",
"Belastbaar inkomen gebouwed","Belastbaar inkomen gebouwd","Verwijzing SAT")
VALUES (GeomFromWKB(?, 3857),?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?)

Other layers, also writing to the same database, operate normally.

Does anyone have an idea where to find a solution?


texts: What's relationship between pixel spacing and resolution?

texts:
Let's take as an example the Sentinel-1 SAR satellite for which data is accessible online.
The pixel spacing in the range and in the azimuth are found inside the xml files bundled with the product. But I need the resolution instead.

For a level-1 SLC product, the spacings are given as:

<rangePixelSpacing>2.329562e+00</rangePixelSpacing>
<azimuthPixelSpacing>1.396573e+01</azimuthPixelSpacing>

But, how do I calculate the resolution based on the knowledge of thePixel Spacing?


texts: How can I find LayerGroup for a given object in Leaflet?

texts:
So, I have a LayerGroup with polylines, points, etc. in some Layers.
A certain Polyline has the handler:

Polyline.on('click', lineHelper);

Can I find out in this handler which LayerGroup the Polyline belongs to?


texts: Use value of Field in Calculation in Modelbilder

texts:
In Modelbuilder I have set up a model, where in one step I do a "Calculate Field" like "round((!OID!-1- %MyField% /15)/452,1)+1" . For the calculation i use the value of %MyField%, which actually is an Value of a specific Field in the first Row of a Table. I tried to retrieve is with "Get Field Value" - this works. But how do I ensure, that I have the value retrievet at the moment when I want to "Calculate Field"?
Of cource I make an "Precondition" connection, but if "MyField" is 0 (Zero), the precondition says False, so the Calculate Field is not executed.
To explain more: I iterate through features (A) in feature class. In the iteration model i create an fishnet (new feature class (B)) inside the feature. Then i add a field to the new feature class (B).Now i need to calculate this field using the formula mentioned above. In the formula i need to use a value from a field of the feature (A). Is there a way? (I use Arcgis 10.2)


texts: QGIS3 Creating spatial index in a standalone script as a part of output vector layer

texts:
I am making a buffer analisys in a standalone script. Is it possible to automatically during the analisys add index to the output vector layer (".qix" file)?

I tried to add after the analisys like in QGIS2:

layer_out = QgsVectorLayer(out, out[-7:-4],'ogr')

layer_out.dataProvider().createSpatialIndex()

where out is a path to shapefile, out[-7:-4] the shapefile name.
However in QGIS 3.2 Bonn, I got the error message:

AttributeError: 'QgsVectorLayer' object has no attribute 'createSpatialIndex'

If an automated way of creating index during analisys is not possible, how to create it after it in standalone script?


texts: Cannot make feature layer in python toolbox using parameters

texts:
This should be so simple, I feel, but I cannot get it to work. The issue I'm having right now is making a feature layer in this python toolbox. When I try to create boxLayer I consistently get an error that the input is not supported or does not exist. I would like it to create a feature layer based on the name and number of the box requested in the tool's parameters. There are a couple of other mistakes here, but this one issue is the one I really need solved. Is there a better way to pass the box number parameter to help make a feature layer of that particular box?

class Toolbox(object):
    def __init__(self):
        """Define the toolbox (the name of the toolbox is the name of the
        .pyt file)."""
        self.label = "Fire Department Tools"
        self.alias = "fdept"

        # List of tool classes associated with this toolbox
        self.tools = [BoxBldgCount]

    class BoxBldgCount(object):
        def __init__(self):
        """Define the tool (tool name is the name of the class)."""
        self.label = "Box Zone Building Count"
        self.description = "The user will enter a box number and a
        building     type." + \
        " The tool will add a field and store the building count for
        the
        selected type."
        self.canRunInBackground = False

This is where I'm setting the parameters

    def getParameterInfo(self):
        """Define parameter definitions"""
        param0 = arcpy.Parameter(
        displayName = "Enter the box number",
        name = "inBox",
        datatype = "GPLong",
        parameterType = "Required",
        direction = "Input")

        param1 = arcpy.Parameter(
        displayName = "Select a building type",
        name = "bldgType",
        datatype = "GPString",
        parameterType = "Required",
        direction = "Input")
        param1.filter.type = "ValueList"
        param1.filter.list = ["Single Family","Multi-Family","Commercial","Industrial",\
        "City Property","Storage Sheds","Schools","Church"]

        params = [param0,param1]
        return params

        boxNumber = parameters[0]
        buildingType = parameters[1]

        # Make feature layers from the user input.

And this is where I'm having the error. Why does that input not create the layer by concatenating the box number to the FireBoxMap? For example, FireBoxMap_202 exists in the geodatabase, but when I enter 202 as input I receive an error.

        boxLayer = arcpy.MakeFeatureLayer_management(r"C:\EsriPress
        \GISTPython\Data\City of Oleander.gdb\FireBoxMaps\FireBoxMap_"
        + str(boxNumber), "layer")


texts:
I have a list of river gauging station locations, with the coordinates given in a local system:

L/24(5.5*8.15)
O/13(3.5*0.15)
P/3(0.9*2.1)
M/17(13.5*1.6)

Presumably these coordinates are pointing to a particular sheet and then the location on that sheet.

Has anyone got some more information on using these coordinates, maybe some sort of key to the map sheets? Ideally I want to convert to something usable (maybe decimal degrees).


texts: Iterate over geopandas dataframe returning nearest point on each LineString

texts:
I'd like to iterate through a GeoDataFrame that contains a column of LineStrings and return a list of points; one from each LineString that is closest to a reference point.

#### .  Import data from Overpass API

import overpy

api = overpy.Overpass()

pHigh = api.query("""
(
node[highway=primary](34.978011,-116.697296, 36.747345, -113.976445);
way[highway=primary](34.978011,-116.697296, 36.747345, -113.976445);
rel[highway=primary](34.978011,-116.697296, 36.747345, -113.976445);
);
(._;>;);
out center;
""")

##### . Make a shapefile

import fiona
from fiona.crs import from_epsg
schema = {'geometry': 'LineString', 'properties': {'Name':'str:80'}}
shapeout = "test.shp"
with fiona.open(shapeout, 'w',crs=from_epsg(4326),driver='ESRI
Shapefile', schema=schema) as output:
    for way in pHigh.ways:
        # the shapefile geometry use (lon,lat)
         line = {'type': 'LineString', 'coordinates':[(node.lon,
        node.lat) for node in way.nodes]}
        prop = {'Name': way.tags.get("name", "n/a")}
        output.write({'geometry': line, 'properties':prop})

##### .  Make a geopandas dataframe

import geopandas as gpd

gpanda = gpd.read_file("test.shp")

##### . Point to measure distance from

from shapely.geometry import LineString, Point

lat=36.069173
long= -115.013759

point=Point(lat,long)

##### .  return the coordinates for the nearest point on the linestring
#####    relative to the reference point ('point').  Doesn't Work!

result=[]
for index, row in gpanda.iterrows():
    result += row.geometry.interpolate(row.geometry.project(point))

##### . Seems to work for individual rows...

gpanda.geometry[1:2].interpolate(gpanda.geometry[1:2].project(point))

This may be the wrong approach, I really just want an accurate measure of the nearest highway to many points. This way I would have the nearest point from each highway, which I could use to calculate the great circle distance from the reference point and find the minimum.


texts: Defra LiDAR - Convert DTM to PointCloud?

texts:
Defra has DSM and DTM data for Blackwell in Birmingham, but no Point Cloud data.  I am trying to get LiDAR Point Cloud data to build a golf course model (in TGC2019) for Blackwell, and this uses ChadsTool, and this tool requires Point Cloud format data.
How can I convert DTM to Point Cloud?  Apparently LasTools might be able to do this, and I can download the exes, but I do not know how to convert.  The DTM files seem to have large TIFF files in them, and these  are probably 'the data', if that helps.


texts: Matching shapefile to CSV within QGIS - all strings

texts:
I am trying to link a CSV to a shapefile using the method as described in this tutorial: http://maps.cga.harvard.edu/qgis/wkshop/join_csv.php

It seems to work fine but the problem is that the CSV when imported in QGIS reads all data as strings. I found a suggestion online to save another file as csvt. In the same location (see Changing CSV layer attribute value from text to number in QGIS?) and I tried using “text wrangler” and  “text edit” but it didn't work out.

Can anyone suggest a solution that works on a MAC (QGIS: 2.4.0-Chugiak)?


texts: Numeric issue while loading floating point-valued GeoTIFF into PostGIS

texts:
I am trying to import a floating point GeoTIFF image into PostGIS (2.1), using the following command:

raster2pgsql -s 4326 -I -M -F -C tif_dir/*.tif public.map | psql -U postgres -h localhost -d opengeo

If I omit the -C option, then there everything is OK. However, when -C is used, I get errors about  numeric field overflow and Unable to add constraint: 'nodata_values'.  Skipping.

In either case, select * from raster_columns returns empty.

The nodata value as reported by gdalinfo at the bottom of this question is -3.4028234663852886e+38 (, which is the most negative value of single precision floating point numbers). Is floating point value of this size (32-bit) not supported by PostGIS.

How should I fix this issue here associated with nodata?

raster2pgsql output:

Processing 1/1: tif_dir/prec1_nr.tif
BEGIN
CREATE TABLE

...
NOTICE:  SQL used for failed constraint: ALTER TABLE public.map ADD CONSTRAINT enforce_nodata_values_rast CHECK (_raster_constraint_nodata_values(rast)::numeric(16,10)[] = '{-3.40282346638529e+38}'::numeric(16,10)[])
CONTEXT:  PL/pgSQL function _add_raster_constraint_nodata_values(name,name,name) line 48 at RETURN
PL/pgSQL function addrasterconstraints(name,name,name,text[]) line 98 at assignment
PL/pgSQL function addrasterconstraints(name,name,name,boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean) line 53 at RETURN
NOTICE:  Returned error message: numeric field overflow
CONTEXT:  PL/pgSQL function _add_raster_constraint_nodata_values(name,name,name) line 48 at RETURN
PL/pgSQL function addrasterconstraints(name,name,name,text[]) line 98 at assignment
PL/pgSQL function addrasterconstraints(name,name,name,boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean,boolean) line 53 at RETURN
WARNING:  Unable to add constraint: 'nodata_values'.  Skipping
...
COMMIT
VACUUM

And, gdalinfo gives the following description of the GeoTIFF in question:

gdalinfo prec1_nr.tif
Driver: GTiff/GeoTIFF
Files: prec1_nr.tif
       prec1_nr.tif.ovr
       prec1_nr.tif.ovr.aux.xml
       prec1_nr.tfw
       prec1_nr.tif.aux.xml
Size is 7200, 3600
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-30.000000000000000,30.000000000000000)
Pixel Size = (0.008333333333333,-0.008333333333333)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -30.0000000,  30.0000000) ( 30d 0' 0.00"W, 30d 0' 0.00"N)
Lower Left  ( -30.0000000,   0.0000000) ( 30d 0' 0.00"W,  0d 0' 0.01"N)
Upper Right (  30.0000000,  30.0000000) ( 30d 0' 0.00"E, 30d 0' 0.00"N)
Lower Right (  30.0000000,   0.0000000) ( 30d 0' 0.00"E,  0d 0' 0.01"N)
Center      (   0.0000000,  15.0000000) (  0d 0' 0.01"E, 15d 0' 0.00"N)
Band 1 Block=128x128 Type=Float32, ColorInterp=Gray
  Min=0.000 Max=337.000
  Minimum=0.000, Maximum=337.000, Mean=6.705, StdDev=18.508
  NoData Value=-3.4028234663852886e+38
  Overviews: 3600x1800, 1800x900, 900x450, 450x225
  Metadata:
    STATISTICS_COVARIANCES=342.561389035033
    STATISTICS_MAXIMUM=337
    STATISTICS_MEAN=6.7054569494428
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=18.508414006474


texts: How to relate feature class and table in ModelBuilder?

texts:
I am beginner in ModelBuilder, so I have some trouble with it. I have feature class, which contains forest stands (polygons) and table with data about age of their etages. For example, one forest stand can have the first etage (50 years) and the second etage (100 years) - relation 1:n. How can I join this table and FC ? I have to create tool, which select forest stands according to age (for example, I want to select all stands, which have 100 years etages) but I am not sure how to create it.

Tool in ModelBuilder with "add join". It works uncorrectly. JPRL=shape with stands, vek=table with ages.

Shape "jprl" with stands (polygons).

Table "vek" with ages (column "vek") and etages (column "ET").

ID field is IDPS.


texts: Issues with OSM Plugin in QGIS on Ubuntu 12.04

texts:
(not sure how to wrap the subject into a question! But...)

Here are some issues I've been having with the OSM plugin:

1) When I bring in the .OSM file, close QGIS, then reopen, I get a "Handle Bad Layers" dialogue...see image below:

2) When I open QGIS with the OSM plugin enabled, it seems to cause a crash with Unity Desktop, that is, the menu bar at the top of Ubuntu goes blank, then seems to reload. Disabling the plugin and loading QGIS solves the issue.

Anyone else had this?


texts: What are the appropriate parameters to convert .asc file to geotiff using gdal_translate?

texts:
My starting file of interest is an ASC file with the following gdalinfo

Driver: AAIGrid/Arc/Info ASCII Grid
Files: sd_1_3_110805.asc
Size is 11880, 9720
Coordinate System is `'
Origin = (-117.900000000000010,33.399999999999991)
Pixel Size = (0.000092592592593,-0.000092592592593)
Corner Coordinates:
Upper Left (-117.9000000, 33.4000000)
Lower Left (-117.9000000, 32.5000000)
Upper Right (-116.8000000, 33.4000000)
Lower Right (-116.8000000, 32.5000000)
Center (-117.3500000, 32.9500000)
Band 1 Block=11880x1 Type=Float32, ColorInterp=Undefined
NoData Value=-9999

To translate to a tiff I have tried:

gdal_translate sd_1_3_110805.asc lj.tif

Subsequently I would like to change the projection to UTM flat to ultimately create a 3d model of flat vs geocentric terrain.

Unless I have specified using a geocentric model to osgdem (VirtualPlanetBuilder) to build the 3d terrain model, my model always comes out mirror image (north is south etc...). So I desire to use a flat terrain model... what parameters should I be using to get to a geotiff with the correct projection to start with?


texts: Different results between individual and batch clip of rasters in QGIS?

texts:
When I cut one raster with clip raster by mask layer I get the expected result.

But when applying batch processing with the same vector file, and the same rasters, the rasters aren't clipped and the values are changed to 0 and NaN?

The rasters are landsat 5 images, and the vector has the same datum.

Edit: I've figured out this is the gdal command used when processing individually:

gdalwarp -ot Float32 -q -of GTiff -cutline "some_shape.shp" -co COMPRESS=DEFLATE -co PREDICTOR=1 -co ZLEVEL=6 -wo OPTIMIZE_SIZE=TRUE "original_raster.tif" "clipped_raster.tif"

And this one is used in batch processing:

gdalwarp -ot Byte -q -of GTiff -cutline "some_shape.shp" -co COMPRESS=NONE -wo OPTIMIZE_SIZE=TRUE "original_raster.TIF" "clipped_raster.tif"

Why the differences?


texts: SR EPSG 3035 ETRS/LAEA Europa to EPSG: 6204 Macedonian State Coordinate System python

texts:
I would like change a raster dataset from EPSG: 3035 (ETRS89/LAEA Europa) to EPSG: 6204 (Macedonia State Coordinate System)

My WKT for the first coordinate system is:

PROJCS["ETRS89 / LAEA Europe",
    GEOGCS["ETRS89",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6258"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4258"]],
    PROJECTION["Lambert_Azimuthal_Equal_Area"],
    PARAMETER["latitude_of_center",52],
    PARAMETER["longitude_of_center",10],
    PARAMETER["false_easting",4321000],
    PARAMETER["false_northing",3210000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","3035"]]

and for the Second SR:

PROJCS["unnamed",
    GEOGCS["Bessel 1841",
        DATUM["unknown",
            SPHEROID["bessel",6377397.155,299.1528128],
            TOWGS84[682,-203,480,0,0,0,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",21],
    PARAMETER["scale_factor",0.9999],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1],
    AUTHORITY["epsg","6204"]]

Since I want to reproject my data, and the projections have different datums, how this process would look like using gdalwarp?

According to this post (How to reproject a raster file in QGIS with datum transformation?)
 gdalwarp is the way to go, in the following form:

gdalwarp -s_srs "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" -t_srs "+proj=tmerc +lat_0=0 +lon_0=21 +k=0.9999 +x_0=500000 +y_0=0 +ellps=bessel +towgs84=682,-203,480,0,0,0,0 +units=m +no_defs" input.tif output.tif

Is that a valid approach? I can't see the datum information in the proj4 strings I am providing for both map projections. Is there anything I have to take care of when projecting a raster datasets from a map projection which uses a different datum to my target Spatial Reference?


texts: How to Handle Proxies within Carto's Python Import API framework

texts:
I am attempting to programmatically upload new data sets (tables) to Carto with the Carto Import API Python SDK. Unfortunately, I am doing this behind a proxy which is blocking me from successfully reaching the Carto API endpoint.

I have experimented with the requests and urllib libraries to successfully circumvent the proxy block, but I am not familiar with how to implement this same proxy handling within the framework of Carto's import python module. I have the proxy address and a username that can successfully connect through it.

The following code is what I am using to attempt to connect and upload a sample data set.

from carto.datasets import DatasetManager
from carto.auth import APIKeyAuthClient

USERNAME = '{username_placeholder}'

USR_BASE_URL = "https://{user}.carto.com/".format(user=USERNAME)

auth_client = APIKeyAuthClient(api_key="{api_placeholder}", base_url=USR_BASE_URL)

Local = r'C:\temp\carto\carto_test.zip'

dataset_manager = DatasetManager(auth_client)
dataset = dataset_manager.create(Local)

The error I receive when running this code is as follows:

  carto.exceptions.CartoException: HTTPSConnectionPool(host='{placeholder}.carto.com', port=443): Max retries exceeded with url: /api/v1/imports/ (Caused by NewConnectionError(': Failed to establish a new connection: [WinError 10060] A connection attempt failed because the connected party did not properly respond after a period of time, or established connection failed because connected host has failed to respond',))

Would anyone be able to assist in adding functionality for proxy handling?


texts:
I am trying to have my openlayers map fit the size of the browser and also resize every time the browser resizes. So I can have responsive design.

My structure of code is like so

<script type='text/javascript' src='OpenLayers.js'></script>
<script type='text/javascript'>

var map;
var options = {
  controls: [],
  maxResolution: "auto",
  projection: "EPSG:900913",
  units: 'm'
}

function init(){
 map = new OpenLayers.Map('map_element', options);
 //set layers, styles for vectors layers, set and activate click events etc.
}

</script>

<body onload='init();'>

<div id='map_element' style=" width:800px; height:400px ; z-index:1; position: relative;text-align:center;"> </div>

According to this I tried using

//if browser loads
//this is inside init() because = body.onload > init > resize map
map = new OpenLayers.Map('map_element', options);
map.updateSize();

//if browser resized
    window.onresize = function()
    {map.updateSize();}

but does not work

Also, if I set style=" width:100%; in the map_element div the map disappears from the page.

What is happening? What am I missing?

Thanks in advance


texts: In ArcGIS, using 'Select by Location' and pushing selected features to a new table or feature?

texts: I have a layer with about 100 buildings and another layer that shows the extent of a flood. When I use the 'select by location,' it selects the ~50 buildings that would be underwater if a flood occurred. Is there a way to take those 50 selected polygons and push it to a new/separate layer or table?


texts: Reprojecting non 4326 WMS image to spherical mercator

texts:
In OpenLayers, I'm trying to display a WMS overlay (EPSG:4269) from a server I don't have control over to a spherical mercator projection.

Is there an equivalent to this:
http://blog.sumbera.com/2010/02/17/overlay-wms-on-google-in-openlayers/

for projections other than EPSG:4326 to EPSG:3857?

Thanks!


texts: Spatial join (points in polygon) but also including the "observation count" of each point in the join_count field

texts: I have point layers of bird observations, and I want to symbolize certain areas by the number of bird observations in each polygon. Each point in the point layer has a "observation count" field associated with it. I want to include this number in the join_count field that already includes the number of points in each polygon. Since I will use that field to symbolize the abundance in each area


texts: What can go wrong in UTM to WGS1984 conversion (shifted data)

texts:
I try to convert .shp files to OSM data with shp2osm and ogr2osm but the output data is shifted.

Semi-Cross-Posting Warning: I have already asked a questions on the OSM help site, but this questions is more theoretical.

The .prj file I have contains the following (line-breaks added for readability) and the source data was saved using ArcGIS:

PROJCS["Adindan_UTM_Zone_37N",
GEOGCS["GCS_Adindan",
DATUM["D_Adindan",
SPHEROID["Clarke_1880_RGS",6378249.145,293.465]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]],
PROJECTION["Transverse_Mercator"],
PARAMETER["False_Easting",500000.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",39.0],
PARAMETER["Scale_Factor",0.9996],
PARAMETER["Latitude_Of_Origin",0.0],
UNIT["Meter",1.0]]

I am using ogr2ogr for the coordinate transformation:

ogr2ogr -t_srs WGS84 -s_srs "ESRI::subcities.prj" "subcities_converted.shp" "subcities.shp"

When I finally import the data in the OpenStreetMap editor, the data is shifted by around 176m, 222°.

Does this mean the source data and the .prj file is plain wrong? Or I am doing something wrong on the way?


texts:
Can I change the order of how my coordinates will display on QGIS?

I would like to display Latitude, longitude instead of longitude, Latitude


texts: Create a non spatial index on an attribute in QGIS before joining two layers

texts:
I’m trying to join a .dbf which contains no spatial data to a shapefile containing spatial data. Ideally I’d like to create an attribute index on the join field in the .dbf before I join it to the shapefile.

Does anybody know how I can do this?

Also in the join menu there are two tick box options:

Cache join layer in virtual memory
Create attribute index on join field

Can anybody tell me what these two options mean? Is the attribute index created during the join if I was to tick this or is it done before the join?

I'm using QGIS 2.0.1 on Windows 7.


texts: “No rows were georeferenced” error in CartoDB

texts:
I am trying to geolocate my data by department name and country, but receive an error:

No rows were georeferenced It seems that some of your rows didn't finish successfully. Perhaps these rows contained empty values or perhaps we just didn't know what the values meant. We encourage you to take a look and try again.

Basically, I have an .xls dataset with columns: Homicides - that contain data on homicides rate; Department - the department where this homicide rate was registered; Province - the province to which each department belongs and Country - Argentina. I've been trying to map the location for the homicide rate, but I did not succeed with even one department. What might I be doing wrong?


texts:
I'm trying to count how many floods occur within a given grid cell in ArcMap and I can't figure it out. I've been reading some other posts on here but I haven't been able to find an answer.

I have flood records for several different years in a city. I want to count how many floods occur where. I made a fishnet grid and overlayed that over my city and points, but now I don't know how to get it to count how many points fall within each grid. Does anyone know what tool can do that?


texts: Summing Attributes from a Related Table?

texts:
I have a line feature class representing road segments, with each road segment having a unique ID. This feature class has an attribute of average observed daily traffic counts. I have another table which contains more information about each road segment, but each road segment ID in the table has 8 entries (each of which corresponds to a 3-hour period).

I've set up a Relate between the feature class and the table, and have added a new field to the feature class that I would like to populate with the sum of 8 values from the table that share the same road segment ID (i.e. one attribute, represented by 8 measurements from different times-of-day).

I'd like to do this using the Field Calculator or Model Builder. Python would be preferred over VBScript.


texts: Generating shortest route from A to B (polygons) using Network analyst using ArcMap ModelBuilder

texts: I am trying to make a model that generates the shortest route from polygon A to polygon B with network analyst using ArcMap ModelBuilder. I have managed to get the model to work using point data (point A to point B) but not with polygons. Please advice


texts: How can I modify a VRTRasterBand sub-class,etc from Python?

texts:
I am using GDAL to merge a series of overlapping rasters using Pixel Functions to combine the color bands. As a proof-of-concept, my approach has been to create a VRT, and then to modify it in a text editor:

Add the VRTDerivedRasterBand subclass to the VRTRasterBand tags
Add <PixelFunctionCode> (and related) tags to implement the functions.
Cut&paste the <ComplexSource> tags between bands. Specifically, the new Green (2) band required input from both Red (1) and Green(2) bands. The order here matters, but as long as it is consistent then my custom function can be modified accordingly.

This proof-of-concept had three input rasters, but I wish to work with a variable number.

I also need to automate the entire sequence with one script or program. I'm currently working with Python, but could switch to C/C++.

Can I use the GDAL API to perform these manipulations of the VRT? Can I change the sub-class?

Looks like SetMetadataItem() might be usable to add new XML tags(?) but how can I change the sub-class of an existing VRTRasterBand? Or to remove and replace with a new one?

I can't create a new blank VRT and build it from the ground up, because
I do not know the grid size and transform for the final image.

Or is the best approach to directly modify the VRT's XML using something like XSLT?


texts: Checking if point is contained by any geometry

texts:
I have a table of geometries (mixed MultiPolygons and Polygons), I want to check if a point is contained by any of the geometries in the table. I don't care about how many, nor I want to retrieve any information, all I care about is a boolean. Is there a way to speed up this query, given the circumstances?

SELECT N>0 as on_land FROM(
SELECT COUNT(*) N FROM lands l WHERE
ST_Contains(l.geom,ST_SetSRid(ST_Makepoint(47.3,19.4), 4326))) B

(the union of the 'lands' table is all landmass on earth)


texts: Why Count points in Polygon QGIS2.99 don't work properly?

texts:
I'm using QGIS 2.99 and try to count point in polygon with the function in vector analysis count points in polygon but don't have a correct result.

It must generate a polygon layer with the field count, but it is no generating all the polygons.

I use the same function in QGIS 2.18.16 and it works fine.

Here is the output from both.


texts: Downloading large attachments from GeoNetwork

texts:
I'm not able to download any attachment which is larger than 1GB, I'm getting the following error
<apiError>
runtime_exception
<description>Handler dispatch failed; nested exception is java.lang.OutOfMemoryError: Required array size too large</description>
<message>NestedServletException</message>
</apiError>

Our catalogue should allow downloads for large files (up to 10GB). Is there a possible way to download large attachments. I'm running GeoNetwork with Jetty server, also the same happened when using Tomcat.


texts: Managing a large number of layers in OpenLayers/GeoExt and Geoserver

texts:
I am building a web application to facilitate the handling of approximately 2500+ historic maps. I think it would be too difficult for our users to search for specific maps in an ordinary layer tree.

Is there some way to implement a layersearch in a GeoExt/Geoserver stack?

I'd be grateful for any idea how to deal with that even if it's not an "out-of-the-box"-solution.


texts: Point clustering based on a certain value

texts:
I'm trying to cluster points based on the sum of a attributes of its nearest neighbors.
So here are the points in the data set as they are:

The objective is that each point has an "evalfield" attribute associated with it. I want to be able to cluster them such that the sum of the clustered points ad up to a predetermined value that I decide.(last column).

To illustrate with an example:
Here are a collection of points , with their attribute on the basis for which clustering is too occur written above them.
(0.5,0.6 etc.)

I want to be able to cluster them such that the sum of this value in a cluster adds up to a predetermined value (1.1 in this case)

EDIT

So I followed the steps and pretty much copied out the code as is to run it in qgis.

import geopandas as gpd
from shapely import geometry as geom

gdf = gpd.read_file('filteredptlst.shp')
exploded_gdf = gpd.GeoDataFrame()
exploded_gdf.crs = gdf.crs

print('starting ...')

for index, row in gdf.iterrows():
    if float(row['evalfield']) == 0.01:
        exploded_gdf = exploded_gdf.append(row, ignore_index=True)
    else:
        for i in range(0, int(row['evalfield']*100)):
            entry = row.copy()
            entry['evalfield'] = 0.01
            exploded_gdf = exploded_gdf.append(entry, ignore_index=True)

print('processing done ... saving')

gdf.to_file(driver='ESRI Shapefile', filename="result.shp")

print('saved')

But I keep getting this error:

Python Console
Use iface to access QGIS API interface or Type help(iface) for more info
execfile(u'C:/Users/USER/OneDrive/Documents/Studies etc/adrchitectural design/thesis/Thesisgisdata/Siteselectionalgs/alg1_point_filtering_explosion.py'.encode('mbcs'))
Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "C:/Users/USER/OneDrive/Documents/Studies etc/adrchitectural design/thesis/Thesisgisdata/Siteselectionalgs/alg1_point_filtering_explosion.py", line 24, in <module>
    gdf = gpd.read_file('filteredptlst.shp')
  File "C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\geopandas\io\file.py", line 19, in read_file
    with fiona.open(filename, **kwargs) as f:
  File "C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\fiona\__init__.py", line 170, in open
    raise IOError("no such file or directory: %r" % path)
IOError: no such file or directory: 'filteredptlst.shp'

Ive already tried changing the project environment and putting in a complete file path to the location of the file.

I couldn't se anything wrong with the code so I ran it separately in Visual Studio ,  but that doesn't seem to work either , there are no trackback errors but none of the print messages which I put I ( to check) arent't displayed either.

P.S. I'm really new to both GIS and Python and there must be something really obvious I'm overlooking here.


texts: Installing pgrouting on Ubuntu 16.04 fails

texts:
I tried to install pgrouting on Ubuntu 16.04 but I got following error message:

CREATE EXTENSION pgrouting;

  could not open extension control file "/opt/PostgreSQL/9.6/share/postgresql/extension/pgrouting.control": file or directory not found

I downloaded the EnterprisDB Postgres(9.6) Bundle, installed postgres ant then installed PostGis using the stackbuilder.
I was able to create the PostGis extension and the topology extension.
Some other extensions like pgrouting fail to install.

Does anyone know why intallation of some extensions fail?

I did the same on Win10 a few days ago an everything worked.


texts: Remove unlisted values on Symbology on ArcGIS Pro?

texts:
I have a layer file which I use as reference for symbology for future layers. It is based on unique values, however there is over 100 unique values, and each one have their own icon.

If I apply this symbology to a small layer, its legend or symbology table will have over 100 items, even if the layer only uses 3-4 symbols.

Is there a way to remove all icons that are not been used in that layer?

ArcGIS Pro has a button called Add Unlisted value, and I'm looking for the exact opposite: remove all unlisted values.


texts: Building proj/pyproj projection string for local coordinate system?

texts:
I'm struggling to come up with a proj string for a local county coordinate system using pyproj.  The system is defined as follows:

Linear unit: International foot (ift)
Geodetic datum: North American Datum of 1983 (1991)
Vertical datum: National Geodetic Vertical Datum of 1929 (NGVD29)
System: Central Oregon Coordinate System
Projection: Transverse Mercator
Latitude of grid origin: 43° 00' 00" N
Longitude of central meridian: 121° 17' 00"W
Northing at grid origin: 0.000 ift
Easting at central meridian: 3,300,000.000 ift
Scale factor on central meridian: 1.00016 (exact)

My current string is not converting correctly

p = Proj(r'+proj=tmerc +lat_0=43.00 +lon_0=-121.2833333 +x_0=3,300,000 +y_0=0 +ellps=GRS80 +datum=NAD83 +k_0=1.00016 +to_meter=0.3048 +no_defs', preserve_units=True)


texts: Turn string data (of GPS coordinates) into a GPS tracker = frozen map. QGIS2.14/Py2.7/Win7/LAN

texts:
I am currently involved in a student project which is to turn string GPS coordinates into a GPS tracker. This is my current work: First, the UDP data is shared through a LAN (No internet available) therefore, I get it with an UDP socket. The frame: (DDS coordinates)

"lat+4542000lon-00019590alt+658"

Then I translate it in DD coordinates so that I can show it on the map using this code:

plane=QgsVertexMarker(canvas)
plane.setColor(QColor("blue"))
plane.setIconSize(10)
plane.setIconType(QgsVertexMarker.ICON_CROSS)
plane.setPenWidth(3)
plane.setCenter(QgsPoint(lon,lat))
plane.show()

With 1 mark it works very well. But I receive new frame every second and I'd like to move this mark simultaneously. That's why I tried this:

#I assumed that the first point is already shown
while n<3600:
    data,addr=sock.recvfrom(1024)
    lat,lon,alt=prepare(data)
    plane.hide()
    plane.setCenter(QgsPoint(lon,lat))
    plane.show()
    n+=1

But when I run this, Qgis is frozen until the last mark appears while I expect to see one mark after another and to run QGIS normally.
How can I fix it?
--QGIS 2.14, Python 2.7, Windows 7 32bits, no internet--


texts: Passing string from one Python class to a QTextEdit field of another class (QGIS plugin)

texts:
I've an input mask that feeds a spatialite database.

I have one Python file which does all the work and two Qt Ui-Files (Ui_InputMask, Ui_Items).

In the main form (class mask_n_functions(QDialog, Ui_InputMask)) I can open a new form (class items_mask_n_functions(QDialog, Ui_Items)) which shows me some check boxes. After I check some I want to pass the information to a QTextEdit of the main form class.

Here is the shortened code:

        ....

class items_mask_n_functions(QDialog, Ui_Items):

    def __init__(self, parent):
        QDialog.__init__(self, parent)
        self.setupUi(self)

        self.pass_items.clicked.connect(self.getItems)
        self.abort_items.clicked.connect(self.closeItems)

    def getItems(self):
        getItemList = [['val1','description1'],['val2','description2']]
        ItemList = []
        for l in getItemList:
            checkboxstring = str(l[0])
            checkboxname = checkboxstring.strip()
            checkbox = getattr(self, checkboxname)
            if checkbox.isChecked():
                ItemList.append(l[1])
        ItemListText = '%s' % ', '.join(map(str, ItemList))
        QMessageBox.critical(iface.mainWindow(),'Input...', ItemListText)

    def closeItems(self):
        self.close()

class mask_n_functions(QDialog, Ui_InputMask):

    def __init__(self, parent):
        QDialog.__init__(self, parent)
        self.setupUi(self)

        global now
        now = datetime.datetime.now()

        global username
        username = getpass.getuser()

        ...

       self.choose_items.clicked.connect(self.itemCallMask)

        ...

    def itemCallMask(self):
        items_mask = items_mask_n_functions(self)
        items_mask.show()

        ...

I googled: "python how to pass variables between classes", but I don't know how to implement this for my plugin structure.

Ui_InputMask contains the QTextEdit I want to fill with ItemListText of the class items_mask_n_functions.

Within one class I set Text to a QTextEdit like this: self.qtextedit_tofill.setText('some text').


texts:
I'd like to change the name of a workbench I published on FME Server. I'm afraid I didn't find a way to do this.
As I want to keep the history of commits I cannot add the workbench as an new file to the repository.

FME Server Version 2019.2.1


texts: Deriving Values from Attribute Table for Value Table in Python Toolbox of ArcPy?

texts:
I would like to derive the values in the value table dynamically based on the given input shapefile. I managed to derive the values but I'm not able to delete or move any of those values if I use the provided buttons in the toolbox.

Assuming, an attribute table (which is provided as parameter 0) has 10 features and in field "purpose" are the values "purpose1","purpose2" and "purpose3" (each of those values occurs several times but I only want the unique values). The value table should (after providing the shapefile) be automatically filled with those three values from above. If another SHP is provided (lets assume with the values "purpose_a", "purpose_b", "purpose_c", than the former values shall be replaced by the latter.

The values in the value table shall be derived based on a field in the provided input shapefile.

What I did so far:

class ExportWells(object):
def __init__(self):
    """Define the tool (tool name is the name of the class)."""
    self.label = "some_tool"
    self.description = "some_tool"
    self.canRunInBackground = False

def getParameterInfo(self):
    """Define parameter definitions"""
    global sl_purpose
    #list with all possible values of purpose
    sl_purpose = [[u"purpose1"],[u"purpose2"],[u"purpose3"],[u"purpose_a"],[u"purpose_b"],[u"purpose_c"]]

    input_shp = arcpy.Parameter(
        displayName="Input Shapefile",
        name="in_shp",
        datatype="GPFeatureLayer",
        parameterType="Required",
        direction="Input")
    input_shp.filter.list = ["POINT"]

    purpose = arcpy.Parameter(
        displayName = "purpose",
        name = "sl_purpose",
        datatype = "GPValueTable",
        parameterType = "Optional",
        direction = "Input")
        purpose.columns = [["GPString", "purpose"]]
        purpose.values = sl_purpose
    params = [input_shp,purpose]
    return params

def updateParameters(self, parameters):
    """Modify the values and properties of parameters before internal
    validation is performed.  This method is called whenever a parameter
    has been changed."""
    global sl_purpose
    if parameters[0].altered:
        get_selection_values = [[shp_values[0]] for shp_value in arcpy.da.SearchCursor(parameters[0].valueAsText,["purpose_field"])]
        sl_purpose_list = list(set([purpose_value[0] for purpose_value in get_selection_values]))
        sl_purpose_modified = [[sl_purpose_item]for sl_purpose_item in sl_purpose_list if sl_purpose_item]
        parameters[1].values = [profile_status_val for profile_status_val in parameters[1].values if profile_status_val in sl_profile_status_modified]
        #1. attempt - assign derived list (without any global variables)
        #parameters[1].values = sl_profile_status_modified
        #2. attempt - assign those values from global sl_purpose which are in the derived list
        #parameters[1].values = [profile_status_val for profile_status_val in sl_purpose if profile_status_val in sl_profile_status_modified]
    return

Unfortunately, the result is not satisfying at all. My first attempt, to simply assign the list of values to the value table works in so far that the values are written in the value table. Unfortunately, the functionality of the buttons (move up and down, delete etc.) is gone. The same behaviour occurs with my second approach. This has another disadvantage, that I have to provide a list with all possible values in the getParametersInfo (or alternatively in the init which doesn't work neither).

The third approach works a bit better (at least, the functionality of the buttons is given) but it has the same drawback, that I have to provide a list with all possible values (I tried to provide an empty list, but that didn't work). The problem here is, that it doesn't update the value list properly in the case of changing input shapefiles.

Does anybody have an idea how to cope with that?

I'm using ArcGIS 10.4.1


texts: Table Select does not work in Model Builder, when using parameter in SQL

texts:
I'm trying to create a model that uses values from a single table, in order to perform a different calculation for each raster in a file.
Each raster is represented in the table by a row that contains the raster unique name (string), and 2 other fields with the numeric values for the calculation (slope & intercept).
I use "Iterate rasters" to go through the files, and "Table Select" where I select the relevant row for the parameters, by using %Name% (from the iterator) = name (the raster name in the table).
Doesn't work... the model gets stuck in the "Table Select" stage.
Does anyone have any idea why?


texts: QGIS: create lines to connect points which are within a distance to each other

texts:
I have a point layer and I'm Looking for a tool wich connects every point of the layer with each other. Ideally, I can give a distance and the number of Connections for each point. The following example is from the ArcGIS tool

In the example each point connects to 9 other points wich are within distance X Meters.

Is there a tool like that in QGIS or a Plugin. Or perhaps it can be done with geometry Generator somehow?


texts:
Anyone have information about where Print Composer is hiding.
I have followed the instructions on the QGIS website which says it is under Projects or the button in the control panel. Neither option is there. I also don't have a File button which it says is another option.
I use QGIS 3.16 version.


texts: How to edit a given version of geodatabase?

texts:
I have a versioned geodatabase. Any of that versions are used only by one user. So then user called "XXX" has a version called "XXX". So then now I want to make a changes on given version, but when I create a code like that:

IWorkspaceEdit edit = null;
    try {
        edit = new IWorkspaceEditProxy(findVersion(getVersionedWorkspace(), "myVersion"));
        edit.startEditing(true);
        edit.startEditOperation();
        IFeatureWorkspace featureWorkspace = new IFeatureWorkspaceProxy(getWorkspace());
        ITable iTable = featureWorkspace.openTable("MyTable");
        IRow iRow = iTable.getRow(56);
        StringBuilder builder = new StringBuilder();
        iRow.setValue(21, "test2");
        iRow.store();
        edit.stopEditOperation();
        edit.stopEditing(true);
    } catch (IOException e) {
        e.printStackTrace();
    }

I get an exception:

  AutomationException: 0x80041354 - Syntax error in 'Microsoft VBScript compilation error'

When I call this method again a get an information like I should use an edit session

  Objects in this object class cannot be updated outside of an edit session.

but how you can see on above code, I start and stop editing and edit operation.

Anyone can give me any hint for that. I also using an IFeatureWorkspace.update method, but I get an information about problem with edit session. I'm not sure what I'm doing wrong, because I can access to database and get rows, but I can't edit them. I am using a JAVA as programming language.

When I change my code into:

IWorkspace iWorkspace = ArcGISServerHelper.getInstance().getWorkspace();
IWorkspaceEdit edit = null;
try {

    IVersionedWorkspace iVersionedWorkspace = new IVersionedWorkspaceProxy(iWorkspace);
    iWorkspace = new IWorkspaceProxy(iVersionedWorkspace.findVersion("myVersion"));

    edit = new IWorkspaceEditProxy(iVersionedWorkspace.findVersion("myVersion"));
    IFeatureWorkspace featureWorkspace = new IFeatureWorkspaceProxy(iWorkspace);
    ITable iTable = featureWorkspace.openTable("MyTable");
    edit.startEditing(true);
    edit.startEditOperation();

I get an exception:

  AutomationException: 0x80004005 - The XML being loaded could not be parsed. [physical name: OWNER.ANNO_TABLE] in 'ESRI GeoDatabase'


texts:
I used r.terraflow.short in GRASS, but it displayed this error:

r.terraflow.short elevation=raster_bulacanElevation \
                  filled=raster_bulacanElevationFlooded1 \
                  direction=raster_bulacanElevationFlowDirection1 \
                  swatershed=raster_bulacanElevationSink1 \
                  accumulation=raster_bulacanElevationAccumulation1 \
                  tci=raster_bulacanElevationTCI1
stats.out: Permission denied
STREAM temporary files in C:/Users/Rhen/AppData/Local/Temp
(THESE INTERMEDIATE STREAMS WILL NOT BE DELETED IN CASE OF
ABNORMAL TERMINATION OF THE PROGRAM. TO SAVE SPACE PLEASE
DELETE THESE FILES MANUALLY!)

What should I do to make r.terraflow.short run? I tried deleting all the contents of the Temp folder but GRASS still displayed the same error.


texts:
I need the name of each layer in an .msd, but I could not find support using arcpy (ArcGIS 10.0).  Is there a method I missed?

In the current workflow, since I am publishing the document, I use the REST service:

descriptionUrl = url + "/ArcGIS/rest/services/" + mapServer + "/MapServer?f=json"
restResponseFp = urllib.urlopen(descriptionUrl)
mapDescription =json.load(restResponseFp)
layerNames = [str(layer["name"]) for layer in mapDescription["layers"]]


texts: How to make multipart feature call outs start at the same point from the annotation text in QGIS3

texts:
I have a map showing multipart features and map text. The text is too large for the features so I want to use call outs but the call out lines extend from the text in what appears to be the shortest path, I would like all the call out lines to extend from the text at the same point. How can this be accomplished?
In the image below I was able to get some features to label the way I intended (i.e. Green Features on the Left) but the Yellow ones I couldn't without moving the text block so far away that it was off the map.


texts: TauDEM from QGIS toolbox "The following layers were not correctly generated"

texts:
I'm having the same problem as described in complete-qgis-uninstallation-to-fix-taudem-issues. However, that question specifically relates to the reinstallation of QGIS to fix the problem. I've tried this and the issue was not fixed. The TauDEM tools run correctly from the command line but not from the QGIS toolbox.

I'm using QGIS 2.18.19 on Windows 10. TauDEM version 5.3.7.


texts:
I have a doubt. In my plugin I use algorithms from processing sextante. For example:

Processing.runAlgorithm("grass:v.to.rast.attribute", None, inputLayer, 0, Indexes, extent, cellSize, -1.0, 0.0001, outPath)

inputlayer = directory of inputlayer
extent = defined extent of the new raster file
cellsize = option for user to choose the cellsize
outPath = output raster directory
and other parameters that doesn't matter here

but my question is in the attribute field (that is defined as Indexes in the example). This parameter is the attribute field to rasterize. So, usually, I used a combobox to define this parameter but in this case I create, internally, a new column with values and I want to use this column. But I don't know how this parameter are defined, as a list, as a stringm...? Does anyone know that? Or how I can access the attribute created by me?


texts: Applying function to all files in directory using PyQGIS?

texts:
I have a problem for apply this function in all files.
This script is about the extraction of the bands from a raster and I want to use it on all the file raster in the folder. It works on the single file but I can not apply to all files from the folder.

Why?

import glob, os, processing
from osgeo import gdal
in_path="C:/.../folder_input/"
out_path ="C:/.../folder_output/"

file_list = os.listdir(path)

file_list= [file for file in os.listdir(path) if os.path.splitext(file)[1] == '.ecw']

for file in file_list:
print (file)

src_ds = gdal.Open(file)

for file in range (1,file.RasterCount+1):
out_ds = gdal.Translate(out_path + '1' + str(file) + 'ecw',src_ds, format='GTiff', bandList=[i])
out_ds=None

AttributeError: 'str' object has no attribute 'RasterCount'


texts: Is there any open-source building squaring tool?

texts:
I have a set of detached building polygons and I can merge them together in order to form a single building polygon. However, when merged output polygon may not preserve the orthogonality of building sides depending on the orientation of the original detached polygons.

In order to make the sides of merged building orthogonal, is there any open-source building squaring algorithm/ tool?


texts: Defining area of influence for origin-destination dataset (where destination is port) using ArcGIS for Desktop?

texts:
Starting from an origin-destination dataset for a given country, where the origins are city (points) and the destination are ports (points), how can I identify the area of influence for each particular port?

Is there a tool in ArcMap which can automate the creation of the polygons for showing which region is more key to one port or another?


texts:
Does anyone know how to mock IGeometry (e.g. using NSubstitute)? I want to use it within an ISpatialFilter, however I always get an exception on the last line of the code beyond.

  COMException: The geometry property of the spatial filter is invalid.
  (-2147216814)

When I use a realworld-instance such as PolygonClass it works, however if I use this and afterwards do some topological operations (intersect) on them I also use real-world code which I want to avoid.

This is what I have so far:

var poly = Substitute.For<IPolygon, IGeometry>();
poly.IsEmpty.Returns(false);
poly.IsClosed.Returns(true);
poly.ExteriorRingCount.Returns(1);
poly.GeometryType.Returns(esriGeometryType.esriGeometryPolygon);
poly.Dimension.Returns(esriGeometryDimension.esriGeometry2Dimension);
poly.Length.Returns(2);
var filter = new SpatialFilter
{
    Geometry = (IGeometry)poly,
    SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects
};


texts: Return Matching Polygons with Leaflet markerclustergroup

texts:
MarkerClusterGroup is a function that clusters points and return their aggregation.

See: example

In my case I have some polygons to which I calculated the centroids. Instead of returning the lines between centroids when hovering on the cluster. I want to return the polygons associated to those centroids. which means that I'll have a highlight on the polygons that have been clustered.

Basically, this is what I did for the time being, it will create the cluster and return the volume. But I want Onhover to highlight the polygons assiciated with the markers(centroids in my case):

var clusters = L.markerClusterGroup({

              iconCreateFunction: function(cluster) {
                var markers = cluster.getAllChildMarkers();
                var volume = 0;
                for (var i = 0; i<markers.length;i++){
                  if (!isNaN(markers[i].value)){
                    volume += markers[i].value
                  }
                }

                return L.divIcon({ html: '<b> Total Volume: ' + volume + '</b>' });
              }
            });

how can I achieve it ?


texts: How to calculate LAI (Leaf Area Index) using Remote Sensing?

texts: I am new in remote sensing, and I would like to ask the community what are some good methods to calculate the Leaf Area Index (LAI) using Remote Sensing Data, for agriculture mostly applications.


texts: Generate centroids for vector polygons based on raster and land-use data using QGIS

texts:
Although QGIS has the ability to generate centroids for polygons, I am not sure how to generate centroids based on land-use development using raster data.

Here's an example of what I'm trying to do:

I have a vector polygon layer, where each polygon represents a distinct zone for analysis.
I have a color ortho raster layer that has the same geometric boundaries of my polygon layer.
I am trying to develop a weighted centroid point layer based on the development density in the raster, yet still related somehow to the polygon it's originally in.
Currently, I do this visually. It is time consuming to do 200+ zones and not as precise if a mathematical model was performing the analysis.
I'm not sure if there is a plugin or a method for measuring the amount of roof space in a raster and using that and their density as a basis to develop centroids for the area

My Question:

How Can I generate centroids for vector polygons based on development density from a raster layer


texts: Sinuosity is a fairly straightforward concept. How can this be calculated in postgis given a series of linestrings that represent the path?


texts: Conversion from .shp o .xml gives invalid data

texts:
I tried to convert an ESRI shape file (jagdgebiete.zip) with ogr2ogr from the GDAL library with this command:

ogr2ogr -f GPX -dsco GPX_USE_EXTENSIONS=YES Jagdgebiete.gpx Jagdgebiete.shp -nlt MULTILINESTRING

But the output is not correct, as it gives e.g. a track-point
<trkpt lat="154616.85700000077486" lon="545566.38300000037998"></trkpt>
which is definitely not in the region in question which is approx. 46-46° N and 13°E.

I think the original unit is in metres from a reference point (see Austria GK M31).

What is the command for ogr2ogr to make the correct transformation?


texts: MBTiles with GeoServer: how to hide white nodata tiles?

texts:
I published MBTiles with GeoServer. When I open it in QGIS I get:

You see that the area around the tiles is white.
Before I generated the MBTiles I tried this solution but I get the same result.

How can I hide these white tiles? And how to find the reason why they appear on the map? Because when I open the MBTiles directly in QGIS I do not see any white tiles.


texts: Computation with dictionary-values in GEE

texts:
I have FeatureCollection in GEE , so for each feature there is data in the properties about total number of pixels, means and more data.

I would like to add new column with computation based on this data in the dictionary, for example,  new column (or key) that will be the result of computation of bt0204/count, or any other computation.

How can I access this key-value data in order to run computations? How can I run computation over values in dictionaries?

My end goal: to be able to calculate the percentage of different pixels values , for example to run bt0204/count*100, and have the result as new key (column).

Edit: I have tried to use this method just to get the number:

var num=dictionary.get('bt0204');

and then print it, but I get null. maybebe because there are many features? I don't know


texts: Can ArcGIS Network Analyst solve multiple vehicle routing problem?

texts: Can network analyst tool solve multiple vehicle routing problem, or is it only solved by using ArcLogistics Route?


texts: How to extract percentage of land cover from a raster according to multiple shapefiles?

texts:
I have a raster map that represents land cover of a region, and different shapefiles (total 120) in the same region. I want to know the percentage of every land cover for each shapefile.

For example, I want to know % of forest, % of cultivated area, % of grassland, and so on, for each shapefile.

I already tried clipping the raster with the shapefiles in Qgis, but is a slow process. I think with R is possible too, but I have not been able to do it yet:

setwd("/home/shapefiles")
shps <- dir(getwd(), "*.shp") #looking for shapefiles

library(maptools)
list.shape<-list()
for (shp in shps){
  list.shape[[shp]]<- readShapeSpatial(shp) #Reading shapefiles
}

#Opening raster
library(raster); library(rgeos); library(rgdal)
r.landcover <- raster('/home/raster')
shape.integrated<-c()
for (k in 1:length(list.shape)){
  shape.integrated<-gUnaryUnion(list.shape[[k]])  # Joining shapefiles
  r.vals[[k]]<- extract(r.landcover, shape.integrated)  # Extract values of the raster according to the shapefile
  percent[[k]]<- lapply(r.vals[[k]], FUN=I NEED AN OPTION THAT CALCULATE PERCENTAGE OR SOMETHING SIMILAR HERE)  # Percentage
}

I can use Qgis, R, GRASS (Open source).


texts: i want to extract aerosol optical depth on land using modis it seems like I need a value from one of these parameters between Corrected_Optical_Depth_Land and Corrected_Optical_Depth_Land_wav2pl  but don't know what one to use can anyone please tell me the difference between these tow or if it better to use another parameter. ?


texts: Removing the gaps between two vector layers and merge them seamlessly

texts:
I have multiple vector files  needed to merge them without any error, but there is gap between one layer and another layer as shown in image 1 respectively. How can I rectify this error?

I merged the layer files but I can see the gap in the layer as shown in Image-2.

PS: All the operations are carried out in QGIS.
How can I resolve this issue?


texts: How to remove the city label in OpenLayers?

texts:
I have a WMS layer like this:

var wms_layer = new OpenLayers.Layer.WMS(
            'OpenLayers WMS',
            'http://vmap0.tiles.osgeo.org/wms/vmap0',
            {layers: 'basic,clabel,ctylabel,statelabel', transparent: true},
            {isBaseLayer: false, opacity: .7}
            );

How do I remove ctylabel or statelabel from this WMS layer?


texts: Estimating road network distances when database queries are expensive

texts:
I have a large number N of points and I'd like to build their distance matrix with respect to a road network.

Unfortunately, I am not able to perform all N^2 queries to build this matrix, so some of these distances have to be approximated, using a limited number of queries together with other information (e.g. crow-flies distances between points).  Are there any well-known methods for doing this?

For example, I could take each point and only look at the distances to some of its nearest neighbors, as well as a few randomly selected points, and then solve some sort of regression model on this, but I would imagine that someone thought of a much more clever scheme at some point.


texts: Is it possible to use qfield on multiple tablets with multiple users, and sync to a single dataset?


texts: Discrepancy between importing same data from CSV and GPS

texts:
I have been importing similar data by using two different methods:

Entering the coordinates manually in excel and save as CSV, then importing in QGIS through delimited text layer. Column A is name, column B is Easting (X) and column C is Northing (Y). WGS is set at 84 for both map and CSV layer.
Importing same data straight from GPS into QGIS

When I use option 1, my coordinate is corrupted and mapped 50 km from the actual location. I would like to find out why as I will create databanks using CSV.


texts: Validating whether shapefile exists using ArcPy?

texts:
I am very new to Python and am trying to complete my final project with not a lot of success.

My code works without functions, the problem is I need it to work with functions. In the main(): function I want the user to enter an .shp file path. Then I have an isValid(inPath): function that I want to call the file path the user entered and validate if the file exists. I have done this in many assignments before with .csv and .gpx files, but I am not having success with .shp files. When I keep everything under the main function, it works no problem.

Here is my code currently:

**# Import Modules:**
import arcpy
from arcpy import env
from arcpy import mapping
import os.path

**# Set workspace**
env.workspace = #this is where I set my workspace
env.overwriteOutput = True

**# Checks if shapefile is valid**

def isValid(inPath):
    if arcpy.Exists(inPath):
        print 'Exists'
    else:
        print ' Does NOT exist'

def main():
    shpFile = raw_input("Enter a Map Document path: ") # User inputs a shapefile
    shpFile = shpFile.replace("\\","/") #replaces backslashes with forward slashes


texts: Inputting Layers in QGIS for PyWPS Processes

texts: I recently downloaded the PyWPS demo files from github here: https://github.com/PyWPS/pywps-4-demo. I'm having some trouble figuring out how to use some of the processes. In particular, processes like "feature_count" require the user to input a "layer". What should be entered in the text box for the process to get my GML layer in QGIS?


texts: Calculating Bike Lane Density for each Census Block Group using ArcGIS Pro

texts:
How do I calculate the Bike Lane Density (existing bicycle facility miles in each block group divided by the area (in square miles) of the respective block group)?
As shown in the image below, I have bike facility miles shapefile containing data in Polylines and I have Block Groups shapefile in polygons and has area information for each block group.


texts: Error occurred decoding the EPSG code EPSG%3A3857 in GeoServer

texts:
My application is throwing the NoSuchAuthorityCodeException error although the layer preview in geoserver is working perfectly, they are not visualizing in the app. I have just checked the logs and found this error

org.geoserver.platform.ServiceException: Error occurred decoding the espg code EPSG%3A3857
    at org.geoserver.wms.map.GetMapKvpRequestReader.read(GetMapKvpRequestReader.java:266)
    at org.geoserver.wms.map.GetMapKvpRequestReader.read(GetMapKvpRequestReader.java:90)
    at org.geoserver.ows.Dispatcher.parseRequestKVP(Dispatcher.java:1545)
    at org.geoserver.ows.Dispatcher.dispatch(Dispatcher.java:692)
    at org.geoserver.ows.Dispatcher.handleRequestInternal(Dispatcher.java:255)

Caused by

Caused by: org.opengis.referencing.NoSuchAuthorityCodeException: No authority was defined for code "EPSG%3A3857". Did you forget "AUTHORITY:NUMBER"?
    at org.geotools.referencing.factory.ManyAuthoritiesFactory.noSuchAuthority(ManyAuthoritiesFactory.java:478)
    at org.geotools.referencing.factory.ManyAuthoritiesFactory.getAuthorityFactory(ManyAuthoritiesFactory.java:457)
    at org.geotools.referencing.factory.ManyAuthoritiesFactory.getCRSAuthorityFactory(ManyAuthoritiesFactory.java:532)
    at org.geotools.referencing.factory.AuthorityFactoryAdapter.createCoordinateReferenceSystem(AuthorityFactoryAdapter.java:778)
    at org.geotools.referencing.factory.ThreadedAuthorityFactory.createCoordinateReferenceSystem(ThreadedAuthorityFactory.java:635)
    at org.geotools.referencing.DefaultAuthorityFactory.createCoordinateReferenceSystem(DefaultAuthorityFactory.java:176)
    at org.geotools.referencing.CRS.decode(CRS.java:517)
    at org.geotools.referencing.CRS.decode(CRS.java:433)
    at org.geoserver.wms.map.GetMapKvpRequestReader.read(GetMapKvpRequestReader.java:258)
    ... 108 more

What could be the problem?


texts:
Is there a way to get QGIS to write out a world file (.wld/.pgw/.jpw) etc rather than an aux.xml.
I see the aux.xml is in general better but it isn't supported by the other software I am using.
I can manually write the .wld based on the aux.xml but that is a bit of a faff.


texts: Invalid layer when dragging tif archive in QGIS

texts:
I'm new in QGIS. I followed this video https://www.youtube.com/watch?v=6adjrQ_VtHQ&t=445s in order to make a bathymetric map, and when I drag the .tif appears this red window

Invalid layer: GDAL provider

  Cannot open GDAL dataset C:/Users/casa/Desktop/Lab/coquimbo.tif:
  `C:/Users/casa/Desktop/Lab/coquimbo.tif' not recognized as a supported
  file format. Raster layer El proveedor no es válido (proveedor: gdal,
  URL: C:/Users/casa/Desktop/Lab/coquimbo.tif)

I am using QGIS 2.18.15


texts: Get coordinates from current view using arcobjects

texts:
I'm wondering if there's a way to get the longitude/latitude coordinates from the center of a data frame? The problem I'm faced with is this: I want to set a bookmark and have ArcMap tell me the coordinates of the center of that bookmark. (Without using a mouse. This all needs to be automated.)
This is akin to hovering one's mouse over the center of a map and reciting the long./lat. numbers that are shown at the bottom of the ArcMap window.

I'm pretty "green" when it comes to VBA programming in ArcMap 10.0. At first I thought I could utilize the "Info" tool, but now I'm guessing that I need to get the use the IEnvelope interface to calculate the center of the map, then grab the coordinates from that. Any ideas?


texts: Simple condition on rasters with different extents, R

texts:
I have multiple raster objects, with the same resolution but slightly different extents. I have use the following function in order to reclassify them and now their values are between 1-6.

What I want to do, is to generate new rasters based on conditions (e.g. rasterC<-rasterA==6 AND rasterB == 6). I tried a simple ifelse, but I get a reasonable error, with regards to their extents. So my question is:

Is it possible to to use conditionals for rasters with different extents (e.g. only for the areas that overlap?)

#from this function, let's say that I get two different rasters (A and  B) with values 1-6.
#How can I create a new one with the areas that have the value 6 on both A and B rasters?

my_reclass <-function(s,p) {
  ifelse(s<0 & p<=0.05,1,ifelse(s<0 & (p>0.05 & p<=0.1),2,ifelse(s<0 & p>0.1,3,ifelse(s>0 & p>0.1,4,
  ifelse(s>0 & (p >0.05 & p<=0.1),5,ifelse(s>0 & p <=0.05,6,NA))))))
}

new_raster<-ifelse(A==6 & B==6)

Apologies for not providing a fully reproducible example, but I do not know how to create from scratch rasters with slightly different extents.


texts:
faced with the situation when ol draws lines if the coordinates are "scattered".
for example.

coords as is:

and sorted coords by date:

bottom image is almost correct except sharp lines up. sometimes the gsm signal may disappear and the coordinates are recorded in the memory of the transmitting device and as soon as a signal appears, the unsent coordinates are transmitted to the server.

how can i fix it? how to make openlayers not draw sharp lines?


texts: GDAL calc difference between two NDVI rasters

texts:
I get a series of NDVI rasters every 2 weeks and need to calculate changes. The rasters overlap, but are different in area. The resolution is also a little different in each raster,

raster 1: -0.34533,  0.34533
raster 2: -0.35001,  0.35001

This can change slightly with every raster.

I am trying to automate the process with python, GDAL. the rasters are between 1gb-20gb in size.

I think i am almost there, but when i try to change the pixel size of the rasters to match, there seems to be a shift, which messes up the calcs.

If i can get the pixel sizes to match, i should then be able to calc the differences.

here is the code for adjusting the pixel size:

WorkingFolder = 'C:/Users/me/Desktop/NDVI/Testing/'

for f in os.listdir(WorkingFolder):
    print (f)
    fName, fExt = os.path.splitext(f)
    if fExt.lower() == '.tif':

        ds = gdal.OpenShared(os.path.join(WorkingFolder,f),gdalconst.GA_Update)
        gt = ds.GetGeoTransform()

        gtLst = list(gt)

        gtLst[1] = 0.36
        gtLst[5] = -0.36

        # apply the changes
        ds.SetGeoTransform(gtLst)

Has anyone done anything similar? or know a better solution for matching big tiffs to subtract? I believe that i am able to match the extents and do all other calculations, but the pixels are driving me crazy


texts: Selecting all values with  by using SelectLayerByAttribute_management with input Parameter

texts:
I am trying to select all values with Null by using the SelectLayerByAttribute_management. But the Fieldname in the where-clause is alterable by using the input-parameter.

That is what i tried:

if inputValue == "Null":
        ap.SelectLayerByAttribute_management(inputFeature,"NEW_SELECTION", inputFieldname +"IS NULL")

Does somebody have an idea what I did wrong?
That's the error message i get:
ExecuteError: ERROR 000358: Invalid expression
Failed to execute (SelectLayerByAttribute).


texts: Conditional statement in field calculator using Python parser ArcGIS

texts:
I have looked at many threads on here asking about coding nested if-then or conditional statements in field calculator in ArcGIS. This seems to be a common problem for beginner programmers. I've attempted to learn some Python and I really wanted to figure this out myself, but am apparently failing.

The Problem:
I have 3 columns (a,b,c) that have numeric values. I want a statement to assign them a value based on their numeric values > 0. For example, If a>0, then 1 and if b>0, then 2, and if c>0 then 3. There may or may not be a situation where 2 columns > 0 and I'm not sure how to deal with that.

I've tried variations of the code below. My errors are incorrect syntax on line 1. Am I using an incorrect function? I can't find any sites listing Python functions that include half of what I see people using here. Am I referencing the wrong columns or the columns incorrectly? I my spacing wrong?

Pre-Logic Script

def Calc( Cockle, Butter, Gaper)
  if Cockle >0:
    return 1
  elif Butter >0:
    return 2
  elif Gaper >0
    return 3
  else 0

ClamType (new field) =

Calc(!Cockle!, !Butter!, !Gaper!)

I tried function Reclass. What am I missing?


texts:
It looks to me that I cannot mix WMS and vector layers' order in OpenLayers 2.

I have an ordered map.layers object and assign z-index values to the layers based on order in the object, for example, map.layers[20] has a higher zIndex than map.layers[19] but lower than map.layers[21].

Yet, if map.layers[20] is a WMS layer and the others are not, the WMS will appear at the bottom on the map, when rendered. The only way to make WMS layer appear on top is to actually setting its zIndex to the highest value on map.layers. Orders between other types of layers are respected and they overlay correctly.

Is it possible to mix WMS and vector layers at all? Or have all the WMS layers to be under all vector layers uncondicionally?

Based on How to define layer order in Openlayers? I manually reordered the whole map.layers object and still it seems to not to work, even if I can check that layer.div.style.zIndex is the theoretically correct for every single layer, and yet that value seems to not be respected when rendering.


texts:
I have some formatting issues with my labels.  I am using a fairly straight forward expression (see below).  It sort of works fine but some formatting is being ignored, namely the ALL CAPS ACP tags. All owner names are lower case with first letter capitalized.  I want to display them as all caps and only SOME labels come out with owners capitalized. I am also using italics  tags when a claim number value exists and that works fine.

Another issue is that although I am  NOT spreading characters in Maplex, some labels but not all have Owner's names spread, but strangely not the PIN numbers and claims numbers all of which are a part of one single string variable.

There are other issues with Stacking Labels in Maplex which doesn't seem to stack properly when using Stacking Separator but let's ignore that for this question.

Is there a way to fix this or do i need to export to annotations and fix manually?

Function FindLabel ([Owner], [MNDMF_Active_Mining_Claims_CLAIM_NUM], [PIN], [PIN_Surface], [PIN_Mining] )
FindLabel = ""

if not isnull([Owner]) then
if trim([Owner]) <> "" then
FindLabel = FindLabel & trim("<ACP>" & [Owner] & "</ACP>")
end if
end if

if not isnull([MNDMF_Active_Mining_Claims_CLAIM_NUM]) then
if trim([MNDMF_Active_Mining_Claims_CLAIM_NUM]) <> "" then
if len(FindLabel) > 0 then
FindLabel = FindLabel & vbnewline
end if
FindLabel = "<ITA>" & "Claim: " & [MNDMF_Active_Mining_Claims_CLAIM_NUM] & "</ITA>"
end if
end if

if not isnull([PIN]) then
if trim([PIN]) <> "" then
if len(FindLabel) > 0 then
FindLabel = FindLabel & vbnewline
end if
FindLabel = FindLabel & trim([PIN])
end if
end if

if not isnull([PIN_Surface]) then
if trim([PIN_Surface]) <> "" then
if len(FindLabel) > 0 then
FindLabel = FindLabel & vbnewline
end if
FindLabel = FindLabel & trim([PIN_Surface])& " S"
end if
end if

if not isnull([PIN_Mining]) then
if trim([PIN_Mining]) <> "" then
if len(FindLabel) > 0 then
FindLabel = FindLabel & vbnewline
end if
FindLabel = FindLabel & trim([PIN_Mining])& " M"
end if
end if

End Function


texts: How to Create Gradient map from GRID file in QGIS-3 Above?

texts:
I have a very important Question. I have a GRD file which can be opened in SURFER or QGIS or GIS.

The problem is that how can I customize and make my grid file perfect?

Here I have two photos which will explain my issue correctly:

this file is the one I'm expecting to create but it's been made using Surfer program. I want to make it this much perfect with QGIS but it only gives me the result below after putting too much time on it:

if there's a better way to do this please explain it with details. QGIS has so many options and I'm not completely familiar with them.

I have created contours from My raster data, then polygonized and categorized them using Graduated type.

the results are shown below:

As you can see the map has become full of square pixels. I need to have it just like Surfer SMOOTHLY as were shown.

Regards.


texts: Open specific page of a PDF using QGIS actions

texts:
Can anyone help me with the syntax to open a .pdf-file at a specific page using a QGIS action? For example, a 50 page pdf has the information of use at page 8, so invoking the action in the QGIS Map Window will open the .pdf at page 8
The qgis Table has 2 columns, one with the location of the .pdf document and one with the page number.
Is this even possible within QGIS?


texts: Number a point layer starting from the furthest north point following a linestring

texts: How can I number a point layer starting from the furthest north point following a linestring clockwise? In the image bellow I have two layers, a line layer and a point layer, and how I would like to number it. Is there a way to do it? My knowledge on virtual layers are very limited to find a solution.


texts: How to Make a New Shapefile Based on Attributes of an Existing Shapefile?

texts:
In my last post I had an issue where I wanted to create multiple shapefiles based on only one attribute field (years). How to do it was answered in that post. I want to make that tool more dynamic so that the user can select any field and any value and make a shapefile for it. Any suggestions on how to use the tool dynamic so that fields and the operator(<,<=,=,>,>=) can also be defined by the user itself.

I tried the following:

import arcpy
fc = arcpy.GetParameterAsText(0)
yrs = arcpy.GetParameterAsText(1)
yrs = yrs.split(';')
subtract = ''
for yr in yrs:
     where = '"arcpy.GetParameterAsText(3)" arcpy.GetParameterAsText(2) ' + str(yr) + subtract
     subtract = ' AND "arcpy.GetParameterAsText(3)" > ' + str(yr)
     filename = str(arcpy.GetParameterAsText(3)) + '.shp'
     arcpy.FeatureClassToFeatureClass_conversion(fc, "arcpy.GetParameterAsText(4)", filename, where)

del yr

I get the following error message:
: ERROR 999999: Error executing function.
An invalid SQL statement was used.
An invalid SQL statement was used.
Failed to execute (FeatureClassToFeatureClass).

Also attaching the screenshot for properties while adding the script.


texts: How do I get the coordinates of a click on vector feature/layer in OpenLayers transformed to another projection system?

texts:
I need to get the coordinates of the click when the user clicks on any place on the base map. I have checked the answer to similar problem here:
How do I get the coordinates of a click on vector feature/layer in OpenLayers?

But, I need the coordinates transformed to the projection system used by my map: that is
EPSG:900913,  Spherical Mercator Projection and passed to a html textbox.

How can I do that?


texts: When I have to use Bing Map Aerial as basemap in ArcMap 10, I get an NVIDIA error, then the main page closes down. How to solve this problem?


texts: Auto-select suitable UTM Zone based on Grid Intersection

texts:
I was wondering if the following approach a) would be considered valid and b) someone already did this in python.

I have a python script where users can choose an area of analysis (anywhere on the world). The tool will select data based on this rectangle (data in WGS1984 lat/lon) and process it. The result is a shapefile (from Shapely) that can be mapped in QGIS or Arcgis (etc.). My current problem is: the users of my script don't have any background in GIS. Therefore, the tool should do as much automatically as is possible. As a further automation step, I would like to output the resulting Shapefile in a UTM-Projected Coordinate System suitable for the chosen area.

users will always choose a pretty small area (local scale), therefore, it is very unlikely that an area stretches significantly across two UTM-Zones
Minimal Distance/Area distortions are important, but only to some degree - the users will perform some calculations that require consistent distances (e.g. GetisOrd Statistic)
 UTM-Grid on Wikimedia
I think I should be able to select the best UTM-Zone by intersecting with the UTM-Grid:

Is there any python package that will detect the best UTM-Zone given an lat/lng point (and perhaps an area extent)? Are there any concerns regarding my approach?


texts: Getting information on geo file projection in Python?

texts:
I have many geo files with different extensions, from different sources (mid/mif, shape, geojson, dat). I'd like to get information on their projection (i've no access to epsg code).

Maybe, I search for the equivalent of ogrinfo linux command in Python. I can't use the following code, because the final program has to be functional on Windows and Linux :

# Python code
args = ['ogrinfo', '-so', '-al', 'destination_data.shp']
subprocess.Popen(args)

The goal is to get srs/crs info as wkt (as the following result) in order to use the http://prj2epsg.org/ api :

GEOGCS["GCS_WGS_1984",
DATUM["WGS_1984",
    SPHEROID["WGS_84",6378137,298.257223563]],
PRIMEM["Greenwich",0],
UNIT["Degree",0.017453292519943295],
AUTHORITY["EPSG","4326"]]


texts:
I'm trying to convert a shapefile into the simple Well Known Text (WKT) format in the example below (from http://en.wikipedia.org/wiki/Well-known_text).

MULTIPOLYGON (((40 40, 20 45, 45 30, 40 40)),
((20 35, 45 20, 30 5, 10 10, 10 30, 20 35),
(30 20, 20 25, 20 15, 30 20)))

I know that I use QGIS to save the shapefile as a CSV, and I can use ogr2ogr to do any numer of conversions. In fact ogr2ogr -f CSV out.wkt source.shp -lco GEOMETRY=AS_WKT gets me as close as I've come so far, but not quite there. Any suggestions?

p.s.  As noted by Mapperz, this thread is very similar to this one. That thread, while it provided the ogr2ogr approach that I noted above, did not solve this particular challenge. The ogr2ogr output I'm getting looks like this (I've truncated the lines). Apparently, I just haven't figured out how to get ogr2ogr to use multipolygon.

WKT,AREA,PERIMETER,PINUPOND_,PINUPOND_I,CODE
"POLYGON ((-120.630531650950061 50.838562484449184, . . .
"POLYGON ((-123.206067372148397 51.038984074378327, . . .


texts: Consolidating Table rows in Postgres/PostGis

texts:
I have a table of weather data as follows:

Table Name: WindGrib

UGRD    -179    -90 2.9
VGRD    -179    -90 5.9
UGRD    -179    -89 1.6
VGRD    -179    -89 6.4
UGRD    -179    -88 -4.2
VGRD    -179    -88 11.3
UGRD    -179    -87 -4.6
VGRD    -179    -87 10.3
UGRD    -179    -86 -12.2
VGRD    -179    -86 13.1

UGRD and VGRD are at the same lon/lats in every case

The source for this data is a wgrib2 extract from a grib file and it looks like this:

"2013-10-21 00:00:00","2013-10-21 00:00:00","UGRD","planetary boundary layer",0,-90,-3

That is all we get ..

The date columns and the "planetary boundary layer" all have the same data.

I need to be able to move the above data into a new table as follows

New table Name: Wind

Geom(POINT,4326)    UGRD    VGRD
-179/-89            2.9     5.9
-179/-88            1.6     6.4

and so on

What PostGres/PostGis Query should I use?