texts:
In my attribute table I have a column named 'BYGN_STAT' which includes information about building status (RA, IG, FA etc.) and dates (YYYYMMDD).

I would like to split this column, and separate building status 'IG' with the following dates. The length of the dates varies from one date to many.

I need to find 'IG' and the following dates before the next to letters. In this case this is 'FA', but it could also be another combination of two letters.

In the figure below I would like to find and copy the values in the red rectangle to a new column. Is this possible?


texts: Convert a .png file to a GeoTIFF while knowing the boundings

texts:
Disclaimer: I've already found this link Convert PNG to GeoTiff using GDAL - but it was of no help.

Actually I've got a .png file of a hand-drawn map (I also got it in the .svg format in case this may help).
I currently know the exact boundings of this file:

   {
     ne: { lat: 48.506161, lng: -1.714582 },
     sw: { lat: 48.477861, lng: -1.742703 }
   }

It's located in the north of France.

So here's the command I tried so far with gdal_translate:

 gdal_translate -of GTiff -a_srs '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'
 -a_ullr 48.506161 -1.714582 48.477861 -1.742703 myfile.png output.tif

My final goal is to serve this file under a tiles system (x/y/z.png) for Leaflet or OpenLayers.
So I used gdal2tiles to convert the output.tif file:

sudo gdal2tiles.py -p mercator output.tif  result/

But my problem is that the generated tiles are geolocated in the north of Madagascar instead of France. Any clue on this?


texts: Mapping paddy rice fields with PMI(Paddy rice mapping index) with Google Earth Engine

texts:
I would like to perform mapping paddy rice with the Google Earth Engine for large scale. I came across with an article which proposed method for mapping paddy rice with special index related to phenological characteristics of paddy rice. Here is link for research paper www.mdpi.com/2072-4292/10/3/447/s1.
According to this article they calculated PMI (Paddy rice Mapping Index) in two crop season, transplantin/seeding and harvest period. Based on this algorithm paddy rice is identified with more than 0.2 thereshold.

Thus, my goal is develop this algorithm in the Google earth engine Platform. I am going to this Multispectral Images such as Landsat and Sentinel-2. Here I wrote some codes but I recieved some error.

Please see my code below,
Map.centerObject(roi);

var spring = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
                          .filterDate('2020-05-20','2020-06-15')
                          .filterBounds(roi)
                          .map(function(image){return image.clip(roi)});

var harvest = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
                          .filterDate('2020-09-10','2020-09-30')
                          .filterBounds(roi)
                          .map(function(image){return image.clip(roi)});

/// PMI Identification

var spring_NIR = spring.select('B5')
var harvest_NIR = harvest.select('B5')

var numerator = harvest_NIR.subtract(spring_NIR);
var denominator = harvest_NIR.add(spring_NIR);
var PMI = numerator.divide(denominator);

var colors = ['white','yellow'];
Map.addLayer(PMI, {palette: colors}, 'Colored PMI');

Code link https://code.earthengine.google.com/2f096e88a9babe58c1e790d64af50642


texts: How to detect corner of shapefile in python

texts:
Is there any way to detect corner of each shapefile in python !? (not BBox)
Or somehow to detect the edges?

I want to detect each shapefile vertices to calculate the differences in latitude for my program.

I've got some idea for that:

Loop over the points going round the polygon and count the number of times the sign of the difference in latitude between point i and i+1 changes.
Calculate consecutive three points angle, then pick four sharpest angles


texts: How does ArcGIS determine which transformation to use when projecting on the fly? Where this would come in would be going from NAD 83 to WGS 84. Multiple transformations could be considered so how is one picked?


texts: Geocoding search bar input into city with Nominatim and Leaflet

texts:
Is it possible to transform a search bar into a geocoder with Nominatim and Leaflet?
Using the search bar I would like to search for a city to view weather information and view it on the map.
In the image below you can see the app I am developing.

In the example I am using geolocation, I would like to implement geocoding to search for a city.
I want to search for the city once to show both its weather conditions and its position on the map at the same time. So I want to turn the search bar into a geocode controller to show all the cities that contain the typed word.


texts: Why are attributes missing in overlay tab QGIS 1.8.0?

texts: I want to add a diagram overlay onto a shapefile layer. When I look at the attribute table, all the correct columns are shown. However, when I open the layer property box, the attribute columns I want to display are not available in the drop down list. Would this be because there are some NULL values in the attribute table? How do I get around it?


texts: Select by attribute how to only select numbers in a string field

texts: Using a SQL experssion in the "Select by Attributes" table, how would I only select records 1, 4, 5, 8 and omit the rest (I only want to select values that are entirely numbers)? What I have right now is SUBSTRING(address,1,1)<='9' AND SUBSTRING(address,1,1)<> ' ', but it won't catch record 2(1234 S Complex).


texts: MultiValue Selection with MultiValue Input

texts:
I made this Python Toolbox as a practice run for a larger project I want to do as I am new to python toolboxes. While this was successful I am not able to the result that I want.

I want to be able to have multiple feature classes as inputs and then allow the users the chance to pick fields they would like to keep from all inputs. Either in check boxes or a drop down list. Much like when you have one input and then you can pick fields using the drop down menu.

I've explored the Field Map datatype to get the desired result but so far unsuccessful.

import arcpy

class Toolbox(object):
    def __init__(self):
        """Deleting Fields as selected from drop down list"""
        self.label = "Delete Toolbox"
        self.alias = ""

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

class DeleteField(object):
    def __init__(self):
        """Delete fields from check list"""
        self.label = "Delete Fields"
        self.description = "Delete fields"
        self.canRunInBackground = False

    def getParameterInfo(self):
        """Define parameter definitions"""

        InputTable= arcpy.Parameter(
            displayName= "Input Table",
            name="InputTable",
            datatype= ["DEFeatureDataset", "DETable"],
            parameterType= "Required",
            direction= "Input",
            multiValue=False)

        Fields = arcpy.Parameter(
            displayName="Fields to delete",
            name="DeleteFields",
            datatype=["Field"],
            parameterType="Required",
            direction="Input",
            multiValue=True)

        Fields.parameterDependencies= [InputTable.name]

        parameters=[InputTable,Fields]

        return parameters

    def isLicensed(self):
        """Set whether tool is licensed to execute."""
        return True

    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."""
        return

    def updateMessages(self, parameters):
        """Modify the messages created by internal validation for each tool
        parameter.  This method is called after internal validation."""
        return

    def execute(self, parameters, messages):
        """The source code of the tool."""

        InputTable  = parameters[0].valueAsText
        DeleteFields   = parameters[1].valueAsText

        des=arcpy.Describe(InputTable)
        fieldlist=des.fields
        DeleteList=[]

        for field in fieldlist:
            if not field.required:
                if field.name in DeleteFields:
                    DeleteList.append(field.name)
        for f in DeleteList:
            arcpy.DeleteField_management(InputTable,f)
            return InputTable


texts: Find out which raster has highest values

texts:
I have 3 rasters for water use and 1 raster for water supply.

I want to know which one has the highest values which will mean that it has the greatest influence on available water supply.

What tools can I use for this in arcmap 10.0 and how do I do this?


texts: Dealing with geometry that should be geography (bad coordinates)

texts:
I'm working on a project with a few publicly available datasets. One of them is from the USGS, it's their prediction of peak ground acceleration in 50 years. Here's a link, and here's what it looks like loaded into QGIS:

The problem is that the "coordinates" of San Diego, California are (roughly) -1,700,000 by 3,000,000.

The real data I'm working with is stored in SQL Geography, in a SQL Server, with SRID 4328 matched up real-world latitude and longitude.

My question is, especially since I've seen a lot of this type of who-knows-geometry going on, how do I mesh up geography with weird geometry? I'm using OGR2OGR to convert from shapefiles, in this case.


texts: Importing GeoJSON into PostGIS with Python

texts:
I am trying to import GeoJSON file with Python to my PostgreSQL db.

Here is full code:

with open('C:/xampp/htdocs/app/data/'+a) as file:
gj=geojson.load(file)
for feature in gj['features']:
    geom=feature['geometry']
    kc_broj=feature['properties']['BROJ_KC']
    ko_broj=feature['properties']['KO']
    ko_naziv=feature['properties']['K.O.']
    jls='RUGVICA1'
    id=str(ko_broj)+'_'+kc_broj
    cur.execute('INSERT INTO dkp (id,geom,ko_broj,kc_broj,ko_naziv,jls) VALUES (%s, ST_GeomFromText(ST_AsText(ST_GeomFromGeoJSON(%s))), %s, %s, %s, %s)', (id, geom ,ko_broj, kc_broj, ko_naziv, jls,))

And I get error :

psycopg2.ProgrammingError: can't adapt type 'Polygon'

Does anyone know how to fix this?


texts: Writing custom BatchProject script using ArcPy?

texts:
My GPS collects data in WGS1984 and I have decided to start working in UTM instead of state plane.  Needless to say, I would edit shapefiles, not paying attention to what coordinate system I was using and relied on the transformations.  I have realized the error of my ways and want to correct my workflow.

I would like to take my original files from my gps and do a batch project to UTM so I can correct the data.  The directory structure of my gps software is such that each project has it's own folder with subdirectories of corrected/uncorrected data and the shapefiles.  I would like to be able to load up the tool with the various files from different locations and save the UTM file to the same directory as the WGS1984 file, just with UTM at the end of the file name.

What I have done is made a copy of the BatchProject script from ArcToolbox to modify, but I have little to no experience with Python.  If this is not the way to approach this problem, please let me know.  If I should be using just the regular project tool and put that inside of a loop of some sort, that's fine too.


texts:
I have a REST KML service that I am trying to import into a PostgreSQL table using Python 3. I haven't ever dealt directly with a KML service before, only JSON, and I'm having some trouble finding a way to consume this format. With JSON I use the response library similar to this:

import requests
import responses
response = requests.get('http://url_or_ip')
data = response.json()

I can then iterate through the data variable and get what I need, but I haven't found anything similar for KML. Is there a similar method to iterate through KML?

For reference, the service I'm trying to consume is located here: http://206.74.144.42/eitms/roadconditions/


texts: Water color significantly different from different relative orbits?

texts:
I've been examining satellite imagery (Sentinel-2) and noticed that water can look significantly different in images from different relative orbits (seems to happen consistently for the examples below).  Does anyone know why this is (this is during summer and not frozen/ice)?  I looked at the angles and it appeared there's reasonable overlap in angles (Zenith?), although the brighter one is a larger angle. Both Sentinel-2A, and around the same time of day, and L2A product.
Obviously water is highly reflective -- is it related to that?  Trying to make sure there isn't something about the data quality I need to know about.


texts: QGIS Plugin with DockWidget and MainWindow

texts:
For my QGIS Plugin, that I started with the QGIS Plugin Builder, I built a QDockWidget and populated it with some functions.

Now I want to add a QMainWindow, since I need more space for further options, but would like to keep my QDockWidget, in case that one only needs the "simple" Version of the plugin.

I now have build a simple QMainWindow with Qt and want to add it to my plugin. And this is the point where I don't know how to continue.


texts: Why is a projection conversion off by factor of SQRT(2)?

texts:
I am attempting to draw a 3 km ring around downtown Minneapolis. Let downtown be defined as Government Plaza (-93.2661, 44.9765). I am using this ring to calculate census tracts within a distance of a point, so I am using the Census Bureau's projection; the census shapes are given in lat/lon coordinates, so I have to convert back and forth to calculate shape intersects on one hand, and measure distance on the other.

In Python, using shapely and pyproj:

import shapely.geometry as shpgeo, pyproj

lon, lat = (-93.2661, 44.9765)

proj = pyproj.Proj(proj="utm", ellps="GRS80", datum="NAD83")

# project into UTM to calculate distance in meters
point = shpgeo.Point(*proj(lon, lat))

# get a radius 3000 meter circle
bufferShape = point.buffer(3000)

# convert from utm back to lon/lat
coords = shp.__geo_interface__['coordinates']
# this works if shp.__geo_interface__['type'] == Polygon
newCoords = [[proj(*point, inverse=inv) for point in linring] for linring in coords]
bufferShape = shpgeo.shape({'type': shpType, 'coordinates': tuple(newCoord)})

Now, let us test this shape, to see if we got what we wanted:

print(bufferShape.centroid)
print(bufferShape.bounds)

POINT (-93.26610002878017 44.97650391330887)
(-93.29303779998378, 44.95738308641549, -93.23916202385492, 44.99562317086678)

The centroid is just about dead on through the conversion to and from UTM; at least accurate enough for my purposes. The bufferShape should be a circle around the starting point; the max_x coord and the centroid's y coord should be the north-most point of the circle. That point is (44.9956, -93.2661); which according to distance checks with third party software is 2.12 km away; not 3 km.

How did this happen? The number 2.12 looked suspicous to me, so I changd this line:

 bufferShape = point.buffer(3000 * 1.41421)

 POINT (-93.26610005756029 44.97650782661781)
 (-93.30419574627818, 44.94946643426645, -93.22800390139928, 45.003546080298214)

The new north-most point is (-93.2661, 45.0035). Multiplying the distance by the square root of 2 in UTM, I get the correct coordinates in lon/lat.

Why is this happening?


texts: Unable to load shapefile as feature class into geodatabase file, with only 3 of the fields

texts:
I'm trying to automate importing census boundary shapefiles to feature class in a Geodatabase File, but somehow could not get it to work with this code:

### ======================== load boundary map ===============================

outLocation = "C:/python/RA506Final/RA506FinalPro.gdb"
outFeatureClass = "RamBoundaryMap"
inputTable = "tl_2010_27123_tract10.shp"

if arcpy.Exists(inputTable):
    # List of field headers to convert
    inFields = ['FID', 'Shape', 'GEOID10']
    fm= arcpy.FieldMappings()
    fieldMap= arcpy.FieldMap()
    for field in inFields:
        fieldMap.addInputField(inputTable, field)
        fm.addFieldMap(fieldMap)
    #creates a fieldmap with name of input field
        vars()[field] = arcpy.FieldMap()
        fm.addFieldMap(vars()[field])
        #fm.AddFieldMap(vars()[field])

    #convert it now
    arcpy.FeatureClassToFeatureClass_conversion(inputTable, outLocation, outFeatureClass,field_mapping=fm)

    # Confirm Ramsey/Hennepin Boundary map shapefiles imported into Geodatabase file
    x = arcpy.ListFeatureClasses("*", "polygon")
    for i in x:
        print i
else:
    print inputTable + " does not exist"

print arcpy.GetMessage

Here is the error message:

  FieldMappings: Error in adding a field map to field mappings <function GetMessage at 0x0C4E8330>


texts: Feature math/ computation in Google Earth Engine

texts:
I'd like to find a way in GEE to do column math. i.e. to calculate percentage of a feature. GEE have band math, but not feature/featureCollection math. for example:
var feat = ee.FeatureCollection([
  ee.Feature(null, {Id: 0, count: 0 }),
  ee.Feature(null, {Id: 1, count: 10 }),
  ee.Feature(null, {Id: 2, count: 20 }),
  ee.Feature(null, {Id: 3, count: 20 }),
]);
// add a column/property as the percentage of count
var feat = ee.FeatureCollection([
  ee.Feature(null, {Id: 0, count: 0 , percent : 0}),// 0/(10+20+20) = 0
  ee.Feature(null, {Id: 1, count: 10 , percent: 0.1}), // 10/(10+20+20) = 0.2
  ee.Feature(null, {Id: 2, count: 20, percent: 0.4 }),
  ee.Feature(null, {Id: 3, count: 20, percent:0.4 }),
]);


texts: Custom marker and properties from geoJSON in Leaflet

texts:
I am having some issues with two things in a Leaflet map 1) custom marker 2) pulling attributes from a geoJSON file and putting it in the pop-up.

I should qualify things, I am really new to Leaflet and web development. I have looked at a number of different examples, but I cannot seem to make it work within my code. Here is a section of the code, which is filtering based on an attribute in the geoJSON, but I would also like to give it a custom icon and attributes in the pop-up.

 var soffParkingMarkerPopup = 'parking';
 var soffParking = L.geoJson(locations, {filter:soffParkingFilter
 }).addTo(map);
 onEachFeature: function soffParkingFilter(feature,layer) {
   if (feature.properties.type === "parking") return true
 };
 soffParking.bindPopup(soffParkingMarkerPopup).addTo(map);


texts: Why does the className option of a Leaflet popup change only the outside of the popup?

texts:
I'm styling a popup using the className option in leaflet.

Here's the entire code:

<!DOCTYPE html>

    <meta name="viewport" content="width=device-width, initial-scale=1">
    <meta charset="UTF-8">
    <meta name="author" content="eddiearni">
    <title>example</title>
    <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.7.3/leaflet.css"/>

and for some reason the styling shows up behind the popup, which is still styled according to leaflet's default CSS:

Shouldn't the .popup CSS  predominate in the cascade?

Why is there a default styled popup over the background?

Obviously the default styles are still applying here, I'd love to learn how to override them, in order to create individually styled popups for different map features.


texts:
I have a CSV file which contains lat long coordinates which represent a car doing laps around a racetrack. At the end of each lap, there is a row in the CSV file which marks the end of the lap and does not follow the same format as the rest of the rows within the file. An example of such row would be "# Lap 1". I am trying to read through the file, pull out the lat long values, and create a polyline feature class which will display each of the laps. I have been able to create a single polyline but have not been able to create multiple polylines. I am trying to write the code so that when it recognizes the "#" at the beginning of the line, it will skip reading the line and start a new polyline.

I keep receiving the error "RuntimeError: Point: Input value is not numeric".

I'm fairly new to Python and I'm not sure where to go from here. I would really appreciate any help or input.

Here is my code as it stands currently:

    import csv
    import arcpy

    #Function that creates polyline
    def addPolyline(fc, array, sr):
        polyline = arcpy.Polyline(array, sr)
        with arcpy.da.InsertCursor(fc, ("SHAPE@",)) as cursor:
            cursor.insertRow((polyline,))
        array.removeAll()

    arcpy.env.overwriteOutput = True
    arcpy.env.workspace = r"C:\GEOG 485\Lesson4"
    gpsTrack = open(arcpy.env.workspace + "\WakefieldParkRaceway_20160421.csv", "r")
    lineFC = "laplines.shp"
    sr = arcpy.SpatialReference(4326)
    arcpy.CreateFeatureclass_management(arcpy.env.workspace, lineFC, "POLYLINE", spatial_reference = sr)

    spatialRef = arcpy.Describe(lineFC).spatialReference

    csvReader = csv.reader(gpsTrack)
    header = next(csvReader)
    latIndex = header.index("Latitude")
    lonIndex = header.index("Longitude")
    timeIndex = header.index("Time")

    #Create an empty array object
    vertexArray = arcpy.Array()

    #Loop through the line in the file and get each coordinate
    for row in csvReader:
        if row[timeIndex].startswith("#"):
            if vertexArray.count > 0:
                addPolyline(lineFC, vertexArray, spatialRef)

        lat = row[latIndex]
        lon = row[lonIndex]
        vertex = arcpy.Point(lon, lat)
        vertexArray.add(vertex)

    #Add polyline
    addPolyline(lineFC, vertexArray, spatialRef)


texts: Best Landsat band for forest classification?

texts: Using Qgis and the Semi-Automatic Plugin for Classification, which are the best Landsat bands to use to classify trees species (in central Argentina)? In Landsat 5, 7 and 8?


texts: Generating Table or chart of point values per polygon from time series - Google Earth Engine

texts:
I'm trying to generate a table with my GNDVI points value from a time series.
I have a feature collection of my polygons, another feature collection of my random points in the polygons. I have a image collection as well and I want a table with the points value in every polygon for every date in the time series. I joined the two feature collections and I don't know how to proceed from here.
How do I do that?
my code:
var fc = ee.FeatureCollection(user_account+folder+shpfile_name);
var S2 = ee.ImageCollection(image_collection)
        .filterDate(first_date, last_date)
        .filterBounds(fc)
        .filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 18);

//-----------------------------------------------------------------------------------
print(S2)
print(fc)           // CRS must be WGS84 (for the shapefile)

//-----------------------------------------------------------------------------------
//Image Expression for calculating GNDVI----------------------------------------------
var band_1 = 'B3'
var band_2 = 'B8'

var addGNDVI = function(image) {
  var GNDVI = image.expression(
  '(B8-B3)/(B8+B3)',
  {
    'B3': image.select('B3'),
    'B8': image.select('B8'),

  }).rename('GNDVI');
  //NDVI index
  return image.addBands(GNDVI.add(ee.Image(0.0)));
  //Add 0.0 value to GNDVI (currently not in use)
};
var t = ee.FeatureCollection.randomPoints(fc, 100)
print(t)
Map.addLayer(t)

// Calculate GNDVIfor images in the Image collection
var S2_withGNDVI = ee.ImageCollection(S2).map(addGNDVI);
print(S2_withGNDVI)

// Show image of mean GNDVI (over image collection) on GEE.
var scene = ee.Image(S2_withGNDVI.mean());

// Prepare palette for image display
var palette = ['FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',
          '74A901', '66A000', '529400', '3E8601', '207401', '056201',
          '004C00', '023B01', '012E01', '011D01', '011301'];
          //palette for GNDVI image on GEE.

// Add GNDVI image and shapefile on GEE platform (min;max are range values)
Map.addLayer(scene.select('GNDVI'), {min: min_v, max: max_v, palette: palette}, 'GNDVI', 1, 0.85);
Map.addLayer(fc.style({color: 'e02406', fillColor: '00000000'}));

//Export GNDVI Time Series into a Scatter Chart and a CSV format
var TimeSeries = Chart.image.seriesByRegion(
S2_withGNDVI, fc, ee.Reducer.mean(), 'GNDVI', 3, 'system:time_start', fc_ID)
    .setChartType('ScatterChart')
    .setOptions({
      title: 'GNDVI Time Series For Shift 1-2 Argentina',
      vAxis: {title: 'GNDVI'},
      hAxis: {title: 'Date'},
      lineWidth: 1,
      pointSize: 4,
      series: {
}});
// Display.
print(TimeSeries);

// Define a spatial filter as geometries that intersect.
var spatialFilter = ee.Filter.intersects({
  leftField: '.geo',
  rightField: '.geo',
  maxError: 1
    });
    print(spatialFilter)

// Define a save all join.
var saveAllJoin = ee.Join.saveAll({
  matchesKey: 'points',
});
print(saveAllJoin)

// Apply the join.
var intersectJoined = saveAllJoin.apply(fc, t, spatialFilter);
print(intersectJoined)


texts:
I am trying use the micromap package in R, but something is not working. This is the code:
library(micromap)
library(rgdal)
library(ggplot2)

setwd ('C:/Users/Jesus/Dropbox/linked maps/Chapter4')
mexico<-readShapePoly("ESTADOS.shp")
gdp<-read.csv('mexico.csv')
gdp$CVE_ENT<-as.character(gdp$CVE_ENT)
statePolys <- create_map_table(mexico, 'CVE_ENT')

lmplot(stat.data=gdp,
       map.data=statePolys,
       panel.types=c('labels','dot','map'),
       panel.data=list('estado','percapita', NA),
       ord.by='percapita',
       grouping=8, median.row=FALSE,
       map.link=c('CVE_ENT','ID'))

The code seems to work but produces the the following image with incomplete information:

As you can see in the image some states are not ploted. Any idea on why my code is not working?
The link to the files is here: https://www.dropbox.com/s/idvi1nlxq3uml3h/mexico.rar
Thanks,


texts: GEE, scatter plot. ERROR: array length mismatch

texts:
In Google Earth Engine i am trying to find scatter between elevation and occurrence but i am getting common number of array vectors
'''
var data=ee.Image("JRC/GSW1_3/GlobalSurfaceWater").clipToBoundsAndScale(geometry);
var dem=ee.Image("USGS/SRTMGL1_003").clipToBoundsAndScale(geometry);
var mask = data.select(['occurrence']).gt(0);
var data = data.updateMask(mask);
var dem = dem.updateMask(mask);
Map.addLayer(dem,'','dem');
Map.addLayer(data,'','data');
Map.addLayer(mask,'','Mask');
var pixelValsY = data.reduceRegion(
    {reducer: ee.Reducer.toList(), geometry: geometry, scale: 30});

var pixelValsX = dem.reduceRegion(
    {reducer: ee.Reducer.toList(), geometry: geometry, scale: 30});

var x = ee.List(pixelValsX.get('elevation'));
var y = ee.List(pixelValsY.get('occurrence'));

var chart = ui.Chart.array.values({array: y, axis: 0, xLabels: x});
print(chart);
'''

Erroe msg is
'Error generating chart: xLabels must be the same length as the 1-D vectors from the array.'
please help me out


texts: How to compute Straight Skeletons using Python?

texts:
Is there a Python package which provides an implementation of a Straight Skeleton algorithm?

I'm aware that the open source (C++) project CGALcontains an implementation but it seems like cgal-bindings doesn't include this CGAL package.

In any case, I would prefer a pure Python implementation which I could modify/extend to suit my needs.

While an implementation that can handle polygons with holes would be preferable, it's not strictly necessary.


texts: Populating Arcgis.com Large Thumbnail Field

texts:
I'm consuming content from arcgis.com to build a custom map gallery.  I see that the response objects have a field called "largeThumbnail", but in my case all the values are null.  Anybody know how to configure content so the largeThumbnail property is populated?

Thanks in advance for any suggestions.


texts: Simplify Polygon - no loss allowed, looking for a concave hull - ArcMap

texts:
I am looking for input on a computational efficient way (avaliable in ArcMap/arcpy) to create concave hull bounding boxes that simplifies the geometry without removing any of the original polygon. Basically an "OUTSIDE" option on the simplify tool.

Constraints:
Cannot remove any of the original polygon.
Has to be computationally efficient.
Does not rely on external downloads or scripts (very restricted online access, no administrator priveledges).
preferably fills holes of "sufficiently" small size.

Background: In my organisation, we have several hundred datasets of varying complexity an size. We use a spatial query tool that lets users know what datasets are present (intersects) a selected area. For this, we use a simplified coverage polygon for each dataset. This polygon can have a maximum of 50000 vertices, which is quite easily surpased on national datasets, and must be a single object (multipart polygon)

What I would like is something like this convex hull polygon, but that follows the edges into the convex areas.
The grey area is a dissolved grid selection, part of my current experimentation. simplyfing that polygon does make a nice simplified version, but includes a processing chain of "select by location" --> Dissolve --> simplify - I feel this can be a quite cumbersome operation on large datasets with thousands of polygons.

So, before I jump into the rabithole fully, is there a way to produce a concave hull efficiently within ArcMap / ArcPy?


texts: FME Search Envelope error when using Geopackage?

texts:
I try to use FME Search Envelope to clip my data smaller. Unfortunately it interrupt all after few seconds.

I am using GeoPackage, is it fatal problem to use FME:s search envelope?

2019-03-14 14:43:17|   0.3|  0.0|INFORM|FME Configuration: Start freeing memory when process usage exceeds 95.48 GB of virtual memory
2019-03-14 14:43:17|   0.3|  0.0|INFORM|FME Configuration: Stop freeing memory when process usage is below 71.61 GB of virtual memory
2019-03-14 14:43:17|   0.3|  0.1|INFORM|Creating writer for format: NULL (Nothing)
2019-03-14 14:43:17|   0.3|  0.0|INFORM|Trying to find a DYNAMIC plugin for writer named `NULL'
2019-03-14 14:43:17|   0.3|  0.0|INFORM|Loaded module 'NULL' from file 'C:\Program Files\FME\plugins/NULL.dll'
2019-03-14 14:43:17|   0.3|  0.0|INFORM|FME API version of module 'NULL' matches current internal version (3.8 20180604)
2019-03-14 14:43:17|   0.4|  0.0|INFORM|Creating reader for format:
2019-03-14 14:43:17|   0.4|  0.0|INFORM|MULTI_READER(MULTI_READER): Will fail with first member reader failure
2019-03-14 14:43:17|   0.4|  0.0|INFORM|MULTI_READER(MULTI_READER): Adding GEOPACKAGE Reader with keyword GEOPACKAGE_1
2019-03-14 14:43:17|   0.4|  0.0|INFORM|Using Multi Reader with keyword `MULTI_READER' to read multiple datasets
2019-03-14 14:43:17|   0.5|  0.1|INFORM|Creating reader for format: OGC GeoPackage
2019-03-14 14:43:17|   0.5|  0.1|INFORM|Trying to find a DYNAMIC plugin for reader named `GEOPACKAGE'
2019-03-14 14:43:17|   0.5|  0.0|INFORM|FME Configuration: Source coordinate system for reader GEOPACKAGE_1[GEOPACKAGE] set to `EPSG:3067'
2019-03-14 14:43:17|   0.9|  0.3|INFORM|Coordinate System `EPSG:3067' parameters: CS_NAME=`ETRF89.TM-35/Fin' DESC_NM=`ETRF89 / ETRS-TM35FIN' DT_NAME=`ETRS89/01' GROUP=`EUROPE' MAP_SCL=`1' MAX_LAT=`70.13' MAX_LNG=`31.59' MIN_LAT=`59.3' MIN_LNG=`19.09' PARM1=`27' PROJ=`TM' QUAD=`1' SCL_RED=`0.9996' SOURCE=`EPSG, V6.11.2, 3067 [National Land Survey of Finland; http:/]' UNIT=`Meter' X_OFF=`500000' ZERO_X=`0.001' ZERO_Y=`0.001'
2019-03-14 14:43:17|   0.9|  0.0|INFORM|Loaded module 'GEOPACKAGE' from file 'C:\Program Files\FME\plugins/gdal/GDAL.dll'
2019-03-14 14:43:17|   0.9|  0.0|INFORM|FME API version of module 'GEOPACKAGE' matches current internal version (3.8 20180604)
2019-03-14 14:43:17|   0.9|  0.0|INFORM|FME Configuration: Using FME Reprojection Engine
2019-03-14 14:43:17|   0.9|  0.0|INFORM|Reprojected search envelope from (550000, 6750000) - (560000, 6760000) in coordinate system `EPSG:3067' to (550000, 6750000) - (560000, 6760000) in coordinate system `EPSG:3067'
2019-03-14 14:43:17|   0.9|  0.0|INFORM|GEOPACKAGE reader: Directive '_ENFORCE_FOREIGN_KEY_INTEGRITY' has value 'YES'
2019-03-14 14:43:17|   0.9|  0.0|INFORM|GEOPACKAGE reader: Opening file 'C:\cygwin64\data\mtk\mtk_muokattu.gpkg'
2019-03-14 14:43:17|   0.9|  0.0|INFORM|FME API version of module 'NULL' matches current internal version (3.8 20180604)
2019-03-14 14:43:17|   0.9|  0.0|INFORM|FME API version of module 'NULL' matches current internal version (3.8 20180604)
2019-03-14 14:43:17|   0.9|  0.0|INFORM|Loaded module 'LogCount_func' from file 'C:\Program Files\FME\plugins/LogCount_func.dll'
2019-03-14 14:43:17|   0.9|  0.0|INFORM|FME API version of module 'LogCount_func' matches current internal version (3.8 20180604)
2019-03-14 14:43:17|   0.9|  0.0|INFORM|Loaded module 'InspectorFactory' from file 'C:\Program Files\FME\plugins/InspectorFactory.dll'
2019-03-14 14:43:17|   0.9|  0.0|INFORM|FME API version of module 'InspectorFactory' matches current internal version (3.8 20180604)
Error running translation.


texts: Undefined symbol error when creating postgis extension

texts:
When trying to create postgis extension, I'm getting an error, no references in google on it.

maptest=# CREATE EXTENSION postgis;
ERROR:  could not load library "/usr/pgsql-9.3/lib/rtpostgis-2.1.so": /usr/pgsql-9.3/lib/rtpostgis-2.1.so: undefined symbol: _ZNK11xercesc_3_013XMLAttDefList14isSerializableEv

RHEL 6 64bit
PostgreSQL 9.3.5
GDAL 1.11.0
postgis-2.1.0

Any ideas?

UPD.also, after some looking around:

[root@ror6ax3 ~]# ldd -d /usr/pgsql-9.3/lib/rtpostgis-2.1.so
        linux-vdso.so.1 =>  (0x00007fffd3bff000)
        libgeos_c.so.1 => /usr/lib64/libgeos_c.so.1 (0x00007f5275679000)
        libz.so.1 => /lib64/libz.so.1 (0x00007f5275463000)
        libm.so.6 => /lib64/libm.so.6 (0x00007f52751de000)
        libgcc_s.so.1 => /lib64/libgcc_s.so.1 (0x00007f5274fc8000)
        libc.so.6 => /lib64/libc.so.6 (0x00007f5274c34000)
        libgeos-3.3.2.so => /usr/lib64/libgeos-3.3.2.so (0x00007f52748a2000)
        libstdc++.so.6 => /usr/lib64/libstdc++.so.6 (0x00007f527459c000)
        /lib64/ld-linux-x86-64.so.2 (0x00007f5276381000)
undefined symbol: _ZNK11xercesc_3_013XMLAttDefList14isSerializableEv    (/usr/pgsql-9.3/lib/rtpostgis-2.1.so)
undefined symbol: _ZN11xercesc_3_013XMLAttDefList9serializeERNS_16XSerializeEngineE     (/usr/pgsql-9.3/lib/r                                                                                tpostgis-2.1.so)
undefined symbol: _ZNK11xercesc_3_013XMLAttDefList12getProtoTypeEv      (/usr/pgsql-9.3/lib/rtpostgis-2.1.so)
undefined symbol: H5T_NATIVE_SHORT_g    (/usr/pgsql-9.3/lib/rtpostgis-2.1.so)

The behavior is same when building another versions from sources.


texts:
I would like to check out data from an SDE database but for performance sake, i prefer to delete all features outside of my survey area.

i've tried to copy features in area from sde to a file GDB then checkout , but the checkin could not be perfect, because i can only synchronize new added features to sde database. so i cannot synchronize modified or deleted features.

i've tried to checkout directly from sde Database then delete all features outside area using Bounding box query , but it take much time once i start arcpad after checkout.

i'm using Arcpy to checkin and checkout.

Any help please.


texts: How to update a Vector layer with WFS protocol, after updating the filter?

texts:
I am using a Vector Layer with WFS protocol to show only a subset of the data. The filter which defines the subset, changes dynamically from the code. How do I refresh the layer, to show only the data that satisfies the new filter?

layer.refresh(); does not work.


texts: How do dissolve polygons from shapefile using open source tools?

texts:
In my current project I'm visualizing data on a map representing some values for each administrative areas of a country. I'm given a shapefile of this country, which I'm converting to GeoJSON for use with PolyMaps. Now I've got a requirement to cluster multiple areas into bigger ones, since there is no data for some of the smaller parts, but aggregated data for big ones can be acquired. Those clusters are predefined (e.g. A+B+C=D), so I figured it would be easier to produce a clustered shapefile and convert that to GeoJSON afterwards.

But how do I create such a cluster, preferably automated and with open-source tools? Any other sensible way is of course welcome too.


texts: How to get the slope from an elevation raster in QGIS?

texts:
I need to extract the slope from an elevation raster in order to do a species distribution model.
It should be quite straightforward, but I can't seem to get it to work.
I am running QGIS 3.10 on a Windows10.
I have tried using the QGIS and GDAL functions for "slope". Both give similar results, shown in the pictures below. The first is of the original topo file, and the second is the output where the entire area is fuzzy and uniform. There isn't even a difference between the sea and the land, and I can confirm, this area is NOT flat! I have tried using different topo files, but it has given me similar results where everything is uniform.

Here is the other raster file I have been using, with its output as well.
I have tried to change the symbology, but it remains the same; that isn't the problem...

I have also tried using the GRASS7 function for "r.slope.aspect" but it gives me an error message (shown below) that I have not been able to figure out... However, even if I get the GRASS7 modules working, I don't think this will solve my problem, as the same thing is happening with both the QGIS and GDAL functions.

I am quite new to GIS and am struggling through every hurdle.


texts: Set Minimum Zoom for View in OpenLayers 3?

texts:
How can I set a minimum zoom level for a view in OpenLayers 3?

map = new ol.Map({
    layers: [imageLayer, vectorLayer],
    renderer: 'canvas',
    target: 'map',
    view: new ol.View2D({
        projection: pixelProjection,
        center: ol.extent.getCenter(pixelProjection.getExtent()),
        extent: pixelProjection.getExtent(),
        zoom: 3.5,
        maxZoom: 4, // there is option for maxZoom, but no minZoom?
    })
});


texts: Reading data from file geodatabase using OleDbConnection?

texts:
I'm able to connect to a file geodatabase using an OleDbDataAdapter and read the attribute table just fine - from within visual studio.

However, when I put the code inside a windows service it crashes (no error, just kills the service).

The line it crashes on is commented in the code below (when executing the reader).

    Dim fgdbPath as String = "c:\esridata\myfgdb.gdb"
    Dim tableName as String = "TableXX"
    Dim conn as OleDbConnection = New OleDb.OleDbConnection("Provider=ESRI.GeoDB.OleDB.1;Data Source=" & fgdbPath & ";Extended Properties=WorkspaceType=esriDataSourcesGDB.FileGDBWorkspaceFactory.1;Geometry=WKB")
    Dim reader As OleDbDataReader = Nothing
    Try
        conn.Open()

        Dim cmd As New OleDbCommand("Select top 2 ID, CITY  FROM " & tableName, conn)
        reader = cmd.ExecuteReader() ' crashes here!

I also tried using an OleDbDataAdapter.Fill(DataSet, tableName) but it would crash on this line as well.

The service will run forever as long as it doesn't hit the above line of code.  The try/catch block doesn't catch the error either.

If I can't get this to work, I suppose I could use arcobjects and a cursor but that's the last resort.

I'm using ArcGIS 10, .net framework 3.5, vs 2008


texts:
What are "Lazy Geometries"? How are they different from other types of geometries?

I'm an experienced Python/Django developer but am very new to the world of GIS, and understand that this is a newbie question. I've been unable to find the answer on google (if anyone can direct me towards a good tutorial site, I'd be grateful).

I encountered the term when reading up on GeoDjango (and it's use of the GEOS library).

Thanks in advance.


texts: QGIS Modeler if/else condition, custom script for Modeler

texts:
I want to make a simple if/else function inside my Model and I wanted to know if it's possible or not.

I want this thing to do that if an Input for 'Segment Layer' is selected that it
doesn't have to execute the 'Segmentation(meanshift)' and uses 'Segment Layer' instead as an Input for the next Process.

the Inputs for the 'Segmentation(meanshift)'are as follows

Update after trying to implement Joseph's answer:
I used the following code and got this error

Update 2, after replacing the parameters I get an error in line 28 this is were 'else:' stands

##Example=name
##BGREN=optional raster
##Segment_Layer=optional vector
##output=output vector

if BGREN is not None:
    processing.runalg('otb:segmentationmeanshift', \
    -BGREN\
    -0\
    -70\
    -70\
    -0.1\
    -100\
    -100\
    -0\
    -0\
    -0\
    -'True'\
    -'False'\
    -1\
    -0.1\
    -'Layer'\
    -'DN'\
    -1024\
    -1\
    -0\
    -output)
else:
    output=Segment_Layer

Update it works :)
after a few try's and with Joseph help this code now works like a charm, the final trick was to replace two optional values with 'None'.

##Example=name
##BGREN=optional raster
##Segment_Layer=optional vector
##output=output vector

if BGREN is not None:
    processing.runalg('otb:segmentationmeanshift', \
    BGREN, \
    0, \
    70, \
    70, \
    0.1, \
    100, \
    100, \
    0, \
    0, \
    None, \
    True, \
    False, \
    1, \
    0.1, \
    'Layer', \
    'DN', \
    1024, \
    1, \
    None, \
    output)
else:
    output=Segment_Layer


texts:
I recently observe a shift in the geographic projection of Sentinel 2 data. Please see the screenshots below for clarification (Granule T33UXP, Band 8 only). The first image was taken on 20160314, the second (shifted) on 20160403. I added a polygon as a reference.
This is unprocessed data as downloaded from the API hub.
S2A_OPER_MSI_L1C_TL_SGS__20160314T171554_A003794_T33UXP_N02.01
S2A_OPER_MSI_L1C_TL_SGS__20160403T152610_A004080_T33UXP_N02.01

Does anyobdy else observe this? Is this a known issue or does it occur just in my particular case?


texts: Data-defined size legend QGIS and Map Unit symbol

texts:
In QGIS, I am making maps of the same project in different scale (global and local). To make my symbol looks alike in all scales, I've choose to use proportional symbol with a size in "map unit".

Therefore I have now a problem with the "Data-defined size legend" as my symbol are just way too enormous in the legend.

How to get this right?


texts: Getting ERROR 000622 when checking for arcpy.DefineProjection_management?

texts:
The purpose of my code is to create a loop that copies shapefiles from a folder to my folder and change its projection 'GCS_WGS_1984'. Each shapefile starts with the year follow by the month and '_rfe.shp' (ex. 200701_rfe.2007) which is why I made the loop the way it is. When the first stage of the loop starts I get a error saying that in my arcpy.DefineProjection_management function, my projection is not valid.

How can I get my projection to be 'GCS_WGS_1984'?

        Traceback (most recent call last):
  File "C:\Users\Elvis\Desktop\Spring\Geog376\lab05\lab05.py", line 32, in <module>
    arcpy.DefineProjection_management(shp, projection)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\management.py", line 8184, in DefineProjection
    raise e
ExecuteError: ERROR 000622: Failed to execute (Define Projection). Parameters are not valid.
ERROR 000628: Cannot set input into parameter coor_system.

    #import modules
    import arcpy
    import os

     # set workspace(This is where all your processing will be done)
     #arcpy.env.workspace = 'S://376/spring14/eherrer2/'

arcpy.env.workspace = 'C://workspace/'

     workspace = 'C://workspace/'
     otherFolder = 'C://folder/'

    backup = []
    for shapefiles in range(1,13):
        arcpy.CopyFeatures_management(outPut+'2007'+str(shapefiles).zfill(2)+'_rfe.shp',inPut+'2007'+str(shapefiles).zfill(2)+'_rfe.shp')
        shp = '2007'+str(shapefiles).zfill(2)+'_rfe.shp'
        describeShp = arcpy.Describe(shp)
        shpSR = describeShp.SpatialReference

        projection = 'GCS_WGS_1984'
        arcpy.DefineProjection_management(shp, projection)
        print 'GCS {0}'.format (shpSR.GCSname)


texts: Converting GML to Geodatabase format by using Productivity Suite software

texts:
I am using ArcGIS Desktop 10.2 version and would like to import MasterMap Integrated Transport Network from Digimap to ArcGIS, So I was following the steps:

Install Productivity Suites v 3.4.1
Import licence
Open OS MasterMap Data Converter.
However, I got an error here about initialization error, which is shown below.

What could the problem be?


texts: Importing shapefile into PostgreSQL/PostGIS?

texts: I'm completely new to PostgreSQL and PostGIS. I finally got both installed and have the PostGIS extension installed. I'm running the PostGIS shapefile import manager. I select my shapefile and click import. It says shapefile import completed but I do not see it in the file in the database. Any thoughts?


texts:
I am using d3.
I normally constructed d3.svg.area().interpolate("basis")
then I want to put a circle on outline of the area. But,
because of  interpolate("basis"). The outline is not
at exact place. So, I can't place circle on
outline of the area.

Do you have any idea about how I can solve the problem ?


texts: How to get display of individual polygons from vector shapefile?

texts:
I am creating a shapefile with adding 20 polgons to it.

        poly = ogr.Geometry(ogr.wkbPolygon)
        for p in range(20):#20 should be replaced by pix_no
            ridx = rind1[p]
            cidx = rind2[p]
            value = array[ridx][cidx]
            Xcoord, Ycoord = self.pixelOffset2coord(ds,cidx,ridx)
            outFeature = ogr.Feature(featureDefn)
            poly.AddGeometry(doPolygon(Xcoord-15, Ycoord-15, paddingX=150, paddingY=150, angle=0, polygonType = 'square') )
            outFeature.SetGeometry(poly)
            outFeature.SetField("VALUE0", int(Xcoord))
            outFeature.SetField("VALUE1", int(Ycoord))
            outLayer.CreateFeature(outFeature)
            outFeature.Destroy()

Now, I want to display the individual polygons in a loop. But I am getting only one polygon in first iteration, two in second and so on. How can I display one at a time?

 for feat in features:
            var = len(QgsProject.instance().mapLayersByName('Hlayer')) != 0
            if var:
                QgsProject.instance().removeMapLayer(Vlout)
                del Vlout

            Vlout = QgsVectorLayer("polygon?crs=epsg:32643", "Hlayer", "memory")
            pr = Vlout.dataProvider()
            pr.addFeatures([feat])

            #val0 = feat["VALUE0"]
            #val1 = feat["VALUE1"]

            #set the polygon border color and make it transparent

            single_symbol_renderer = Vlout.renderer()
            symbol = QgsFillSymbol.createSimple({'color':'255,0,0,0','color_border':'#DC143C','width_border':'0.8'})
            single_symbol_renderer.setSymbol(symbol)

            Vlout.triggerRepaint()
            Vlout.commitChanges()
            Vlout.updateExtents()

            #Add the layer to map registry
            QgsProject.instance().addMapLayer(Vlout)

            layers =  QgsProject.instance().mapLayers()
            map_canvas_layer_list = [l for l in layers.values()]
            self.canvas.setLayers(map_canvas_layer_list)
            self.canvas.setExtent(Vlout.extent())

            self.canvas.refresh()
            self.canvas.show()


texts:
I have a problem with QGIS that I've had both with version 2.2 and now with 2.8. I can't open the projects that I have saved. It is a project of 4 layers: a vector with a region of Madagascar, a river layer, a layer joined to the first one and a layer with GPS points.

When I open the project, go out a window with four columns: Layer Name, Type, Provider, Datasource. I click OK and get this message: "There are still 4 unhandled layer(s), that will be lost if you close now."

So I can't open my project.

Does anyone know what the problem is? What should I do?

I have a Mac Osx 10.9.5, 64 bit.


texts: Search directory of shp files for an empty attribute value

texts:
I'm trying to search through a directory of shp files for any shp file with a specific field that is empty. Here is my code as of right now.

import arcpy
from arcpy import env
env.workspace = r'C:\Workspace\Data_State\MS\Administrative_Boundaries\2013_complete'

fc = arcpy.ListFeatureClasses('*.shp')
fields = "MP_TYPE"

for field in fields:
    #get list of values
    values = [r[0] for r
              in arcpy.da.SearchCursor (fc, field)
              if r[0] != None]

    #Check if list is empty
    if not values:
        print field, "is empty"

I keep getting a runtime error 'in_table' is not a table or a featureclass
I've looked at other examples on here and can't get any of them to work.

Ideally:

the script would cycle through the directory
print the shp file with the empty attribute value
then move on  because I'm going to have many directories to search through.

Unfortunately, I can't put the shp files into a gdb for this. I'm using 10.3, so if possible I'd like to use arcpy.da.SearchCursor. I just can't seem to get it to work with the shp files.


texts:
I have some old VB scripts which an I am wanting to turn them into Python.
My old VB code reads
If [LAND_TYPE]="RESVE" Then val = [LAND_NAME]
elseIf [LAND_NAME] = "UCL" Then val = [LAND_NAME]
elseIf [LOTNO] <> "" then val = "Lot " & [LOTNO]
else
end if

LotNo_F = Val

I have attempted to turn this into a Python code, i have used update cursor
fc = CAD
fields = ["LAND_TYPE", "LAND_NAME", "LOTNO", "LotNo_F" ]

# Create update cursor for feature class
with arcpy.da.UpdateCursor(fc, fields) as cursor:

    for row in cursor:

        if row[0] == "RESVE":
            row[3] = row[1]
        elif row[1] == "UCL":
            row[3] = row[1]
        elif row[2] != "":
            row[3] = "Lot" + row[1]
            cursor.updateRow(row)

It all works fine except for my last elif statement, it runs but it does not produce the same output as the VB script.
I am a little stumped to what i am doing wrong?


texts: Create a list/sequence of numbers from two DOY in Google Earth Engine

texts:
I would like to create a list of numbers from a DOY to another DOY every 7 days, dependent on a specific date that I set. This then matches an imageCollection I have, which I then use to create a weekly average.
So far I have this.
var startDate = ee.Date('2015-01-01');

var start = ee.Date(startDate.advance(-2,'month'));
var end = ee.Date(start.advance(1,'year').advance(-1,'day'));

var sdoy = start.getRelative('day', 'year');
var edoy = end.getRelative('day', 'year');

I'm trying to use something like
var DOYList = ee.List.sequence(sdoy, edoy , 7);

But because ee.List.sequence doesn't know to go forward to 365 and back to the edoy, I'm producing a list without any numbers. Is there a way of telling ee.List.sequence to do this? Or another way of creating a list of numbers to use.
My code is here and it's lines 111-118 in the code which are the issue.


texts: Raster calculator doesn't convert values properly

texts:
I have some issues regarding QGIS Raster Calculator in my PC (I have installed two versions-2.18.26 and 3.14.15).
The problem is, I can't further edit rasterized layers which I converted from a vector format. I had set no data values to be 0, though, by using the information icon (i), these values are still appearing not to have data.
A friend of mine (QGIS specialist) told me to use the ("raster">=0)*"raster" function, which should have solved this problem. Nevertheless, this function doesnt change anything in the layer values when I'm executing it through my PC (in neither QGIS version), even though it perfectly works through my friend's PC.


texts:
I have Idea about ArcGIS, and a bit of QGIS and zero programming knowledge, which I am trying for thesedays. I saw GRASS GIS and have started working with shapefiles, which are not being loaded... :( I saw the toolbox and could not find v.in.shape tool. Using file management -> import vector into grass-> v.in.ogr.qgis I cod load the file but could not go farther. I now wanted to access this via python console. Any guide on how should I carry on with my task through which I may be able to use the tools. I can get good help in programming. I tried to get my shapefile read in python console in QGIS but could not.
I am very basic in this so please give easy solutions.

thanks a lot.


texts: Can gdalwarp reproject from ESPG:4326 (WGS84) to UTM?

texts:
I've a GeoTIFF image in ESPG:4326 (WGS84) and when I try with gdalwarp to reproject to UTM, it has problem

gdalwarp -s_srs "EPSG:4326" -t_srs '+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs' -of GTIFF test.tif test1.tif

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

Here is the information of this TIFF file

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 = (-91.773493181429174,29.826496486509463)
Pixel Size = (0.001452475638877,-0.001452475638877)
Metadata:
  AREA_OR_POINT=Point
  TIFFTAG_DATETIME=2015:02:02 14:54:05
  TIFFTAG_IMAGEDESCRIPTION=IDL TIFF file
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=Adobe Photoshop CS6 (Macintosh)
  TIFFTAG_XRESOLUTION=100
  TIFFTAG_YRESOLUTION=100
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( -91.7734932,  29.8264965) ( 91d46'24.58"W, 29d49'35.39"N)
Lower Left  ( -91.7734932,  29.4982370) ( 91d46'24.58"W, 29d29'53.65"N)
Upper Right ( -91.3973020,  29.8264965) ( 91d23'50.29"W, 29d49'35.39"N)
Lower Right ( -91.3973020,  29.4982370) ( 91d23'50.29"W, 29d29'53.65"N)
Center      ( -91.5853976,  29.6623667) ( 91d35' 7.43"W, 29d39'44.52"N)
Band 1 Block=259x10 Type=Byte, ColorInterp=Red
Band 2 Block=259x10 Type=Byte, ColorInterp=Green
Band 3 Block=259x10 Type=Byte, ColorInterp=Blue


texts:
Though there are many answers to this question on the Web, I've found none that are recent or correct.

The prevailing view seems to be that one exports a KMZ from Google Maps which one then edits in Google Earth where it is supposed to be possible to simply drag an element from one folder to another. This doesn't work for me.

There's no way in Google Earth to move an element from one folder to another that I can see.

What is the (currently) correct way to combine layers in Google Maps?


texts:
I have this geotiff file of hillshade which I am trying to render on HTML using Geotiff.js https://geotiffjs.github.io/
but it doesn't seem to render anything.

If I inspect chrome console, it does print console logs without any error but image is just black. Is there a mistake in parameters?

Here is my script
<!DOCTYPE html>
<meta charset="utf-8">
<head>
    <script async src="https://cdn.jsdelivr.net/npm/geotiff"></script>
    <script src="d3-marching-squares.min.js"></script>
    <script src="https://d3js.org/d3.v4.min.js"></script>
    <script src="https://d3js.org/d3-scale-chromatic.v1.min.js"></script>
</head>
<body>
<script>
(async function() {
    let origin = [700, 600];
    let size = [500, 450];
    const tiff = await GeoTIFF.fromUrl("url");
    let image = await tiff.getImage(0);
    const width = image.getWidth();
    const height = image.getHeight();
    const tileWidth = image.getTileWidth();
    const tileHeight = image.getTileHeight();
    const samplesPerPixel = image.getSamplesPerPixel();

    let rasterData = await image.readRasters({window: [origin[0], origin[1],
                                                       origin[0] + size[0], origin[1] + size[1]],
                                                       samples: [0]});

    rasterData = rasterData[0];

    let data = new Array(size[1]);
    for (let j = 0; j<size[1]; j++){
        data[j] = new Array(size[0]);
        for (let i = 0; i<size[0]; i++){
          data[j][i] = rasterData[i + j*size[0]];
        }
    }

    let intervals = [-2, 0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20];
    let bands = rastertools.isobands(data, [0, 1, 0, 0, 0, 1], intervals);

    let colorScale = d3.scaleSequential(d3.interpolateBuPu);

    let canvas = d3.select("body").append("canvas")
    .attr("width", 680)
    .attr("height", 500);

    let context = canvas.node().getContext("2d");

    let path = d3.geoPath()
                 .context(context);

    bands.features.forEach(function(d, i) {
      context.beginPath();
      context.globalAlpha = 0.7;
      path(d);
      context.fill();
      console.log(context)
  });

})();

</script>
</body>


texts:
I am a student at college and recently listened Colorado Public Radio about Meth Possession Arrests. I am interested in obtaining the data if the map creator would like to share it with me for GIS class project about crime in Denver.

Is this possible? Is there a way to leave a question/request for a map author in CartoDB?


texts: How to create a line feature which contains the overlying polygon attributes?

texts:
I have a road network Shapefile which intersects a polygon feature class. I need to
create a excel table which contains the road number, the attribute from the intersecting polygon feature class and the kilometer from which the polygon intersects the line.

Imagine a line 10km long which intersects 3 polygons between the start to the finish. I need a table which tells me the following:

ROAD   From KM    NAME
L242      0      Smith
L242     3.4     Jones
L242     6.2     Andrews

I have ESRI Advanced licenses with network analyst and FME.

Would this be a job for the Linear Referencing toolset?


texts: Leaflet map loads incorrectly when inside a Bootstrap .collapse element

texts:
I'm trying to create a Leaflet map which is hidden when the page loads, and revealed when the user clicks a button.

I'm trying to use the Bootstrap collapse functionality to reveal the map on a button press - sample here.

This works to show/hide the map's div - but the formatting of the map is broken, with a large empty area with no data:

This only occurs when the map's div has the .collapse class - if you remove that class the map loads correctly.

So the question is, how to allow a Leaflet map to load correctly when it isn't visible when the page loads?


texts: How to render features that are just outside the map canvas extent in QGIS 3.4?

texts:
I am trying to generate a cloudy mask, following the recipe in "Qgis Map Design 2nd ed. Edition".  For the most part, I have it working.

But when panning around, the mask changes.  And I think that is due to features just outside map canvas extent not being factored into the size expression calculation. Or the rendering of the the entire layer at that extent.

Here's an example of an area where it should be obscured, but isn't:

And if I pan a bit over, the area is now obscured because those features are now being rendered:

Is there a way to configure the layer to render features that are a specific distance outside the visible extent.  So that features just outside the visible extent with large symbols are drawn on the canvas?


texts: Splitting DXF layers into Shapefile layers in QGIS?

texts:
When I load a DXF into QGIS all the layers are present in the attribute table. Is it possible to split this "entities" layer in QGIS into its component parts and save each component part as an individual shapefile.

I am at present selecting each layer and right clicking and saving as .shp. But is there a fast way because this way is time consuming and boring?


texts: What are my options for creating a tile cache via a WMS server which requires http authentication?

texts: I'd like to create a tile cache and offer the tiles via TMS/WMS. The source is a WMS server that is protected by basic http authentication, but no SSL. I'll also need to reproject the tiles. This is for my own personal use. I'd prefer something that can run on either Windows or Linux, and open source is a huge plus. From some searching it appears MapServer can access WMS servers requiring authentication, but doesn't cache. GeoServer can cache WMS, but can't use authentication. Lastly, TileCache can't reproject. What are my options?


texts: Find sublines, given initial and final offsets, via the command line

texts:
I have a CSV file that looks essentially like this:
id,id_wkt,start,end,wkt
11,11,0,43,"LINESTRING(3.3623676 51.3132714,3.3626663 51.3130068,3.363068 51.3126811,3.3636323 51.3122504,3.3640754 51.3119614,3.3643705 51.3117777)"
12,11,43,217,"LINESTRING(3.3623676 51.3132714,3.3626663 51.3130068,3.363068 51.3126811,3.3636323 51.3122504,3.3640754 51.3119614,3.3643705 51.3117777)"

It consists of records with a WKT linestring. Some of the WKTs appear multiple times and the actual geometry has to be obtained from the linestring using the start and end offsets (length in metres along the line) in their respective columns.
I would like to extract sublinestrings based on these length offsets. Is this possible via a command line tool like OGR?


texts: Custom log message in GeoServer through POST request

texts:
I'm trying to send a string message from my JavaScript client, to GeoServer, and have it show up in the GeoServer logs.

How might I do this?

I turned on 'Verbose' logging, and set the number of characters to log for incoming post requests to 1024, from the global settings page. However, I don't see the message showing up in the log.

My goal is to build a performance testing script, which parses the GeoServer log file, and compares render speed, among other things, for various data configurations.
Here is a code snippet, JavaScript and OpenLayers 5:

map.on('moveend', function(e) {
  d = new Date();
  beginTime = d.getTime();
});

map.on("rendercomplete", function(e){
   if(beginTime > oldBeginTime){
     oldBeginTime = beginTime;
     d = new Date();
     finishTime = d.getTime();
     var elapsedTime = finishTime - beginTime;
     var postData = "Time: " + elapsedTime + " ms";
     var xhttp = new XMLHttpRequest();

     var url = "localhost:8080/geoserver/www/home.html";
     xhttp.open("POST", url, true);
     xhttp.setRequestHeader('Content-type', 'text/plain');
     xhttp.send(postData);
   }
 });


texts: When I use the threshold method to extract water , why are there many areas that are incorrectly classified as water?

texts:
I am using Sentinel-1 radar data. This is a misclassified area:

What caused this, can I improve it? This is my code.
    geometry = /* color: #0b4a8b */ee.Geometry.Polygon(
        [[[73.9753904127242, 54.03780906107912],
          [73.9753904127242, 16.93758362646061],
          [137.6082029127242, 16.93758362646061],
          [137.6082029127242, 54.03780906107912]]], null, false);
dataset = dataset.filter(ee.Filter.eq("ADM0_NAME","China"));
// Load region defined by polygon and add it to the map
//Load Sentinel-1 SAR collection and filter according to data collection type
var S1 = ee.ImageCollection('COPERNICUS/S1_GRD')
  .filterBounds(dataset)
  .filterDate('2016-01-01','2016-12-31')
  .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))

//Add first image to map to get an idea of what a SAR image looks like
//Map.addLayer(S1.first(),{bands: 'VV',min: -18, max: 0}, 'SAR image')
// Filter speckle noise
var filterSpeckles = function(img) {
  var vv = img.select('VV') //select the VV polarization band
  var vv_smoothed = vv.focal_median(100,'circle','meters').rename('VV_Filtered') //Apply a focal median filter
  return img.addBands(vv_smoothed) // Add filtered VV band to original image
}

// Map speckle noise filter across collection. Result is same collection, with smoothed VV band added to each image
S1 = S1.map(filterSpeckles)

//Add speckle filtered image to map to sompare with raw SAR image
//Map.addLayer(S1.first(),{bands: 'VV_Filtered',min: -18, max: 0}, 'Filtered SAR image')
//Here we are using -16. This is only an approximation and will result in some errors. Try adjusting the
var classifyWater = function(img) {
  var vv = img.select('VV_Filtered')
  var water = vv.lt(-16).add(1).rename('Water')  //Identify all pixels below threshold and set them equal to 1. All other pixels set to 0
  water = water.updateMask(water) //Remove all pixels equal to 0
  return img.addBands(water)  //Return image with added classified water band
}

//Map classification across sentinel-1 collection and print to console to inspect
S1 = S1.map(classifyWater).select('Water');
var extracted=S1.median().clip(dataset);
Map.addLayer(extracted,{min: 0, max: 2,alette:['white','black']},'water');
Export.image.toAsset({
  image:extracted,
  description: 'sentinel_Threshold_2016_china',
  scale:30,
  region:geometry,
  maxPixels:1e13,
})


texts: QGIS reprojection from WGS84 to UTM has unexpected 'towgs' in project string which messes up R calculations

texts:
I've been reprojecting some points data to UTM 54S from WGS84 in QGIS. I then open the shapefiles in R and change the units from m to km (because I don't know how to do this in QGIS) - as follows:
boundary <- readOGR(paste(getwd(), "/boundary_1col.shp", sep = ""))
boundary2 <- spTransform(boundary, "+proj=utm +zone=54 +south +ellps=WGS72 +units=km +no_defs")

Then I can use the SpatialPointsDataFrame as an input for a regression model (R package: INLA). Having units in km improves the stability of the process.
This is easy enough when the shapefiles have the projstring "+proj=utm +zone=54 +south +ellps=WGS72 +units=m +no_defs" as revealed by raster::crs(boundary). But today I reprojected one of my shapefiles into UTM in QGIS in the usual way, and loading it into R it came up with the following CRS:
"+proj=utm +zone=54 +south +ellps=WGS72 +towgs84=0,0,1.9,0,0,0.814,-0.38 +units=km +no_defs"

I think the problem was that the shapefile was loaded into QGIS from a .csv file with lat/long coordinates. It was recognized as WGS84 though, and before reprojection the CRS showed as "EPSG:4326 - WGS 84 - Geographic" in Properties. And the UTM reprojection looks innocent: "EPSG:32554 - WGS 72BE / UTM zone 54S - Projected".
So what went wrong....? Something's happened, since my INLA model throws an error when I use the fitted model to make predictions. But how can I correct this? (in QGIS or R, either solution would be gratefully received!)


texts: Using ogr2ogr to convert a DXF to GeoJSON and place in PostgreSQL table

texts:
Is it possible to use either the ogr2ogr command line to both convert a DXF file to GeoJSON type and put into a PostgreSQL table at one time?
e.g.
ogr2ogr -f "GeoJSON" path/to/output.json path/to/input.dxf PG:"dbname=my_database user=postgres" "source_data.json"


texts: Batch visualization of raster (grid) files

texts:
I have about 100 raster (grid) files (in the "ascii grid" format) of drought index. All cover the same area of interest but at different times of year. I need to create images from the grids with exactly same color scheme: index>100 - dark red, index 50-99 - red, index 25-49 orange, index 10-25 yellow, index 0 - 10 white.

I tried to do this in SAGA GIS, but SAGA GIS always automatically generates the color scheme by dividing the full range of values in equal intervals.

Is there any free GIS that you would recommend for doing this type of task? (I have some small experience with bat file or python scripting..)


texts: How should ODC databases be restored and snapshot in Docker?

texts:
I have been trying to make an ODC index database state (products and their corresponding datasets) available (1) as a Docker image - probably kept locally due to its size - or (2) as a Docker volume. I've been restoring databases with compressed SQL files from pg_dump and I want to avoid pushing images that are 10s of GiB in size.
I tried using docker commit to implement option #1, and when running a container based on the resulting image, the size of the /var/lib/postgresql/data directory is the same as in the origin container and datacube system check perceives the connection to the container to be a valid one - so (1) the connection to the origin container and the container from the image from docker commit both appear to be valid to an ODC client - with datacube system check not requesting datacube system init to be run in either case and (2) have the same amount of data in Postgres (same data directory size).
However, for connections to the origin container, datacube product list shows the products and ODC queries in Python return the expected results. In contrast, for connections to the derived container (from docker commit), datacube product list shows no products and ODC queries perceive there to be no products as well.
I also tried option 2 - creating a Docker volume and using the same Postgres image for both - terminating the origin before launching the other container to use the data in the volume. The results are the same.
The ODC client for the origin and later Postgres containers have the same execution environments - same image, same launch settings - but they are in different docker-compose environments (that is, the first docker-compose environment sets up the container with the data, then shuts down, then the second docker-compose environment starts another Postgres container using the same image).
How should this be done? How can databases be restored and then reused elsewhere as Docker containers?


texts:
The Global Positioning System (GPS) is a free service that is owned and operated by the U.S. Government and is always available.

When we buy a GPS device, we don't pay a monthly fee or pay a tax for GPS support. We only pay the price of the device.

I don't understand how it's possible. How does the U.S. make money?


texts: I need an R code to select which polygon in a shapefile has got the minimum distance from a given point.


texts: Python Geoserver gsconfig create resource from Postgis store table

texts:
I've working Workspace and PostGIS store in GeoServer. I would like to use python and gsconfig to create new resource/new layer from existing PostGIS table, with:

computing bbox,
adding abstract,
adding keywords,
setting style,
other attributes.

Is it doable with gsconfig? Could you provide any examples?
The same goes with removing resources (layers).


texts: How to determine drainage area upstream from point?

texts:
I do not know how to come up with the drainage area (km squared) upstream of Lat 77.84, Long -30.96. I have converted my contour shapefile to a raster, identified and filled the sinks and then created a flow direction raster, but after that I am not sure what steps to take.

I can create flow accumulation, but none of this is showing me where to insert my lat and long point to get the drainage area upstream from that point. I have all the hydrology tools, but do not know what to input to get the correct answer. I am a newbie.

I have tried watershed/subwatershed delineation. The results I got were very small areas adjacent to my points. What I need is the whole upstream area that contributes water to that point.


texts:
I'm trying to convert .tif file to .geojson using following cmd:
gdal_polygonize.py flowdirection.tif -f "ESRI Shapefile" flowdirection_shp.shp

After which I'm getting error:
Program 'gdal_polygonize.py' failed to run: Application not foundAt line:1 char:1
+ gdal_polygonize.py flowdirection.tif -f “ESRI Shapefile” flowdirectio ...
+ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~.
At line:1 char:1
+ gdal_polygonize.py flowdirection.tif -f “ESRI Shapefile” flowdirectio ...
+ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    + CategoryInfo          : ResourceUnavailable: (:) [], ApplicationFailedException
    + FullyQualifiedErrorId : NativeCommandFailed

After that I just gave the same cmd but using -q :
gdal_polygonize.py -q flowdirection.tif -f "ESRI Shapefile" flowdirection_shp.shp

Still
Program 'gdal_polygonize.py' failed to run: Application not foundAt line:1 char:1
+ gdal_polygonize.py -q flowdirection.tif -f "ESRI Shapefile" flowdirec ...
+ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~.
At line:1 char:1
+ gdal_polygonize.py -q flowdirection.tif -f "ESRI Shapefile" flowdirec ...
+ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    + CategoryInfo          : ResourceUnavailable: (:) [], ApplicationFailedException
    + FullyQualifiedErrorId : NativeCommandFailed

All other Gdal cmds are working but it seems some problem with gdal_polygonize.
Also I tried gdaltindex , but it gave me some error:
gdaltindex flowdirection_shp.shp doq/*flowdirection.tif
Creating new index file...
ERROR 4: doq/*flowdirection.tif: No such file or directory
ERROR 4: doq/*flowdirection.tif: No such file or directory
Unable to open doq/*flowdirection.tif, skipping.


texts:
I'm looking for an efficient way to filter my 700.000 features by attribute in PyQGIS. As suggested in this post  It looks like using QgsExpression for such a job is a fastest way? Or is there a faster approach? I'm filtering using integers in indexed columns on my database.
Now when it comes to API I've noticed there is a :
for feature in layer.getFeatures(QgsFeatureRequest(QgsExpression.createFieldEqualityExpression('field_name', value))):

do something with feature
Any speed up suggestions?  Is using list comprehension with filtering faster way to go about this problem?


texts: Applying low pass/median filter over DEM in GeoTIFF format using ArcGIS Desktop?

texts:
I want to apply a median/low pass filter (over a few hundred kilometers) over my DEM to reduce the noise. I tried ArcGIS focal statistics tool but no luck. I got a blank GeoTIFF file.

How can I do it and get reasonable results in GeoTIFF?


texts:
I tried all the day to zoom by default on my polygon 'geofield'. I'm using geoserver to store my data, I can display the map but the centre is in the middle of nowhere (I want it to be in the middle of the geofield). I tried something like this:

map.zoomToExtent(geofield.getDataExtent());

but it doesn't work! I also tried with a button like this:

var zoomToExtentControl = new ol.control.ZoomToExtent({
        extent: [-11243808.051695308, 4406397.202710291,-4561377.290892059,6852382.107835932]});

it works when I clicked on it but when I try to replace the coordinates by geofield.geometry.getBounds() the map doesn't display.

Here you can see the full code:

    <!doctype html>

  <link rel="stylesheet" href="http://openlayers.org/en/v3.18.2/css/ol.css" type="text/css">

  <title>OpenLayers 3 example</title>

Weeds map test

Another thing; if you have a wms link for the world map (not satellite view) I'll need it!


texts: SelectLayerByAttribute_management not selecting anything

texts:
I've read through most of the other similar questions and nothing so far has addressed my issue; most are about syntax.

I'm trying to select by attribute from a Feature Layer. But every time, Arc announces that it succeeded and nothing is selected. It works just fine using the same expression to select from a Table View. Any ideas?

expression = arcpy.AddFieldDelimiters(tempFeatureLayer, "BLDG_ID") + " LIKE '{0}%'".format(county)

arcpy.SelectLayerByAttribute_management(tempFeatureLayer, "NEW_SELECTION", expression)
if int(arcpy.GetCount_management(tempFeatureLayer)[0]) > 0:
    arcpy.DeleteRows_management(tempFeatureLayer)

It's an SDE Database feature class, if that makes a difference.

I've also tried rewriting the expression as below with the same result:

expression = """ "BLDG_ID" LIKE '37019%'"""

and

expression = "BLDG_ID LIKE '37019%'"


texts:
I am trying to process some terrestrial LiDAR data so I can utilize forestr CRAN package. The paper for forestr says the vegetation area index was computed independently outside the package. I am just having a hard time figuring out how to calculate the VAI with terrestrial LiDAR data.
I have used CRAN lidR to do a lot of the pre-processing steps for forestr like DTM, normalization, voxelization ext.  I just can't figure out the VAI. If anyone has any experience with TLS calculating do you mind sharing any R code or walking me thru it.

The paper is JW Atkins, G Bohrer, RT Fahey, BH Hardiman, T Morin, A Stovall, N Zimmerman, CM Gough (2018) “ Quantifying forest and canopy structural complexity metrics from terrestrial LiDAR data using the forestr R package” Methods of Ecology and Evolution 9 (10), 2057-2066. DOI: 10.1111/2041-210X.13061
So for TLS he has steps outlined in the supplemental which i think most of it can be done cran package lidR. I just learned how to use R not to long ago. I figured how most of the steps in the paper from LIDR https://jean-romain.github.io/lidRbook/
In the supplemental materials he says "TLS data are imported using the process_tls function. VAI in TLS data must be processed outside of forestr as outlined previously due to inherent specifics of each instrument". process_tls takes an x, y, z¸ VAI format as input with x and y describing the ground position and z describing the height above the ground. VAI then is the VAI for each voxel described by the x, y, z position..
VAI is calculated for each 1 m2 x, z bin. For unsaturated column, VAI is calculated from bin return density, the ratio of LiDAR hits in that x, z bin to the total number of canopy hits in that x-column. For each saturated column, bin return density is adjusted using the values of φ derived from the negative log function that accounts for the probability of canopy elements above a saturation level described in equation S4. A maximum value of VAI (max.vai) is used to derive an observed VAI for each column. It is crucial that this maximum value be a site maximum, not a site mean value of LAI...
This is a small sample of the output I was able to make following all the steps. Is it just math formula to calculate VAI? Or is it more complex?


texts: I am trying to do some annotation in a layer (layer 1) where I draw many single polygons within a specific ROI. To do this I toggle a second layer (layer 2) on and off to compare it to the "basemap" (layer 3 at bottom). Is there a way to set a hotkey to toggle layer 2 while still being in the editing mode of layer 1 (so layer 1 is selected and I want to toggle a non-selected layer).


texts:
I have a huge Geotiff file , and i've tiled it ..
I've created the tiles using gdal_retile.py tool , but i'm getting the following error from Geoserver while linking the layer with the tiles folder : Could not list layers for this store, an error occurred retrieving them: This plugin accepts only File, URL and String pointing to a file any suggestion?


texts:
I am trying to implement 3D barriers from a feature class into a network using a model, but I can't seem to get the resulting barrier to appear on the proper three dimensional plane. The feature class that the barrier is generated from a 3D feature, and displays properly in ArcScene. Yet when I run my model, the barrier that appears in the route layer doesn't retain the 3D information. I need to have the 3D information, as I don't want to apply the barrier to all features just to the features on that specific plane.

I have the barriers set up via a second "Add Location" tool, as a polygon barrier after the Add Location tool that generates the stops. Is there something that I need to set in the tool itself to retain the Z information from the feature class? I do not have "Use Network Location Fields instead of Geometry" selected.

Below a picture of the problem. The green is the feature class, red is the barrier that is created after running the model.

I should also add that I do not need the barrier to be extruded, or represent a three dimensional box or object, just that it needs to have the correct elevation.

I am using ArcGIS 10.1, with the Network Analyst extension, model creating the 3D routes via Model Builder.

Thanks,

eric


texts: Obtaining NDVI data from red-edge imagery?

texts:
I am working with an aerial image containing three spectral bands (blue, green, and red-edge). I can obtain an NDVI image with the imagery using Erdas Imagine, but I do not know which bands Imagine is actually using to calculate the NDVI.

How do I find out exactly which bands Imagine is using to calculate the NDVI?


texts: How to select lines/polygons in QGIS CAD tool?

texts:
I'm intending to create parallel and orthogonal lines by using the CAD tool in QGIS.

For each tool there is a button to select the feature before executing the tool.

But though I fixed all the snapping properties - as far as I think - correctly, the tool will neither snap nor select the focused lines or polygons.

Is there anyone who has to face the same Problem and could help me out of this?

In the screenshot the selection button for the parallel lines tool is indicated by the red arrow.

the missing Picture i mentioned before...


texts: Deleting line vertices inside polygon using ArcPy?

texts:
Is there a way to automatically delete vertices in line features using a polygon file in ArcGIS so that it wouldn't cut up the line but just straightens the line inside the polygon?

The lines are lines and the box represents a polygon that covers them. I want to delete line vertices that are inside the polygon so the line would be straight inside the polygon boundary like on the illustration.

I use ArcGIS 10.1 but have access to other software also if something offers the solution.


texts:
I have been doing the r.walk analysis but fail to generate the result. The movement direction map and walking cost map from stop points is extremely small, limited to a 3x3 or 5x5 square. Does anyone know about this issue?
The DEM is a 90m SRTM DEM.


texts: Unable to select table in "reclassify by table" algorithm in QGIS Processing Modeler

texts:
In the Processing Modeler (Graphical Modeler) for QGIS 3.4.3, the "reclassify by table" tool seems bugged. I cannot select a Reclassification table, even though it's an obligatory field.
Am I doing something wrong here?


texts: Conversion from shapefile to geodatabase feature class in ARCGIS

texts:
I have converted feature class from a geodatabase to shapefile, and after that I have updated some polygons in shapefile. Now I want to get back to the same feature class or a fress/new feature class in the same geodatabase.

How Can I perform this task in ArcGIS ?

The Problem is: I have a defined domains(codded-value-domain) in geodatabase (a dropdown will come for updation of attributes while editing(dropdown contains different land use class names and their codes)), but after converting from feature class to shapefile,in attribute table I only got codes for the landuse classes. After that I have updated some polygons in shapefile with same codes.

How can I get back to the same geodatabase as a feature class which has same domains using either its direct method or some script in ArcPy?


texts:
I am trying to loop through a bunch of mxds in a folder and make some changes, in this case in particular I am trying to replace a text string. This should change all instances of "9/26/2014" with the dynamic date tag. The code runs fine but when I open one of the mxds, the date is not updated. Is there some saving or refreshing that has to be done that i'm not doing currently?

import os
import sys
import string

folderPath = '[folder path goes here]'
sort_list=[]
oldText = "9/26/2014"
newText = "<dyn type='date' format='short'/>"

for filename_1 in os.listdir(folderPath):
    sort_list.append(filename_1) #append to empty list
    list_length = len(sort_list)  #list length for progress
sort_list.sort() #sort

for filename in sort_list:
  fullpath = os.path.join(folderPath, filename)

  if os.path.isfile(fullpath):   #check if its a file
    basename, extension = os.path.splitext(fullpath)  #split into the name and extension
    if extension.lower() == ".mxd":  #if its an ArcGIS document
      mxd = arcpy.mapping.MapDocument(fullpath)
      print mxd
      df = arcpy.mapping.ListDataFrames(mxd,"*")[0]
      print df

      #Find all page layout text elements
      for elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
        if oldText in elm.text:
          print(oldText)
          elmText = elm.text.replace(oldText, newText)
          elm.text = elmText
          print(elmText)
          arcpy.RefreshActiveView()
          mxd.save

          arcpy.RefreshActiveView()


texts: Converting an Esri GDB to GPKG using PyQGIS

texts:
I have been working on putting together a simple script that converts a file geodatabase to a geopackage.  I have a geopackage with empty tables set up so that the data in the gdb can be appended to easily.  I have tried different approaches and want to keep the code simple.  I cannot for the life of me figure out how to make this work though.
Does anyone have any thoughts about what I am doing wrong?  Basically, I am trying to iterate through every featureclass in a GDB and export it/append it to a gpkg.
outputGpkg = r"D:\export.gpkg"
inputGdb = r"D:\import.gdb"
fileEncoding = "UTF-8"

gdbFcs = [ "Feature1", "Feature2", "FeatureN" ] # eventually I want to find a routine to pull fc names from gdb and populate the array

for fc in gdbFcs:
    gdbLyr = QgsVectorLayer(inputGdb + "|layername=" + fc, fc, "ogr")
    print(gdbLyr)
    #if gdbLyr.type() == QgsMapLayer.VectorLayer:
    #    _writer = QgsVectorFileWriter.writeAsVectorFormat(gdbLyr, outputGpkg, fileEncoding, gdbLyr.crs())
    _writer = QgsVectorFileWriter.writeAsVectorFormat(gdbLyr, outputGpkg, fileEncoding, gdbLyr.crs())


texts: Calculate Positional Accuracy of Feature with CE/LE Values

texts:
I would like to report the predicted positional accuracy of features that have been extracted from WorldView imagery. I have predicted horizontal accuracy and vertical accuracy values (CE90/LE90) for each point.

Is there an accepted method for reporting the overall CE/LE for a feature or set of features?

Can you simply use the max CE/LE, or is there a more rigorous method?

The only reference I've been able to find is this paper from ASPRS: http://www.asprs.org/wp-content/uploads/2015/01/ASPRS_Positional_Accuracy_Standards_Edition1_Version100_November2014.pdf

It seems to be focused on sampled accuracy for raster products. It's not clear to me if it's relevant to predicted accuracy of vector features.


texts: SFCGAL compilation: How to specify CGAL path?

texts:
I'm a Linux newbie but nevertheless I want to install SFCGAL on OSGeoLive Virtual Machine (this is Ubuntu 14.04.4. LTS trusty). So first of all I had to get the dependencies, e.g. CGAL. Unfortunately there is only CGAL 4.2 available (at least 4.3 required) "in the packages", therefore I had to compile a newer version from source. But now I have CGAL 4.2 AND 4.6 whereas the latter is located in /home/user/src/CGAL-4.6.3/ and I don't know where is 4.2...
Thus, as expected, when trying to compile SFCGAL with cmake . && make && sudo make install I get the following error:

CMake Error at cmake/Modules/FindCGAL.cmake:32 (message):
CGAL 4.3 is required (found 4.2 in /usr/include/CGAL/version.h)

I thought about specifying the CGAL_DIR manually, as mentioned at https://oslandia.github.io/SFCGAL/installation.html, but it didn't work (maybe just because I don't know how to specify this correctly).

Can somebody please tell me, how to specify the path of "new" CGAL (whether via cmake or as global variable)? Or maybe is it possible to replace the older CGAL version with the newer one?


texts: When generating polylines from a series of points, I used "Points to Line" tool. From this, Arc generates a polyline feature class, which includes the attribute shape_length. What units are these shape lengths in? I am using the Canada_Albers_Equal_Area_Conic projection, which has the units of meters, but the shape lengths do not appear to be in meters. Any tips?


texts: Is there a way to edit the host name of Postgis layer?

texts:
Recently I change my server to a remote host. That's means all my work there I used Postgis layers, the host name at the source is set to 'localhost' and I need to change to 'remotehostname'.
The content of the source appears with the complete connection, host, dbase, login and password and the query with, at the field "source layer" at general tab in the layer properties. I can select copy and paste, but this field can't be edited.
There is a way to edit this source?


texts: setting up php_mapscript in portable GIS

texts:
I'm new here, and I'm also new in webgis, I have download portable
gis, and I love this new version with lot more tools, I have a
question regarding to mapserver that comes with portable GIS v2, how
to set up or add php_mapscript extension to php in xampplite, I have
try to do like this but failed :

Add the necessary files (libmap.dll, gdal17.dll, etc.)  to a new
folder (foldername = A ). I get those files from new version of MS4W.
Add the folder path to the windows path environment.
Modified php.ini to add php_mapscript

When I try to launch apache server it saya something like "couldn't
load msvcr90.dll", I copied the file from windows folder and add it to
folder A, new error says "wrong msvcr version".

Can anyone help me, so I can use php_mapscript extension.

Regards,
Zery


texts: Does charactermarkersymbol work differently in linux?

texts:
I have just started working in GIS. i am trying to project an aeroplane glyph font (character index=111) using charactermarkersymbol but unable to do it in linux.
I am getting the same spherical projection for different character index values.
I am able to create a simplemarkersymbol and all the shapes provided like tetra, cone,diamond are projected.
But, for various font styles all i am getting in charactermarkersymbol is a sphere. below is the code snippet.
  I am working in a linux machine and i'm not sure whether linux supports all these glyph fonts or not. Thanks in advance for your help.

   try {

                IRgbColor color = new RgbColor();
                color.setRed(255);
                color.setGreen(0);
                color.setBlue(0);
                SimpleMarker3DSymbol sym =new SimpleMarker3DSymbol();
                sym.setStyle(esriSimple3DMarkerStyle.esriS3DMSTetra);

                StdFont font = new StdFont();
                font.setName("ESRI Default Marker");
                font.setSize(36);

                CharacterMarkerSymbol characterMarkerSymbol = new CharacterMarkerSymbol();
               characterMarkerSymbol.setAngle(90);

                characterMarkerSymbol.setColor(color);
               characterMarkerSymbol.setFont(font);
               characterMarkerSymbol.setCharacterIndex(111);
                characterMarkerSymbol.setSize(36);

                IMarkerSymbol imSym = characterMarkerSymbol;

                imSym.setColor(color);
                imSym.setSize(500);

                IMarkerElement mkElem = new MarkerElement();
                mkElem.setSymbol(imSym);


texts:
Is there a way to produce a certain amount of random points over different classes in a raster in QGIS?
I have extensions such as GRASS or GDAL, ...


texts: Are there any Amazon Silk demo videos showing tiled mapservice performance?

texts:
After reading about Amazon's Silk browser, it seems like this might work really well serving map tiles that are stored in S3.  Does anyone know of any video from a live demo showing a web mapping app running in EC2 using S3?

Preferably it would be an app with a public url that I could compare performance on a different device.


texts: How open database file *.SQLITEDB in QGIS?

texts: how open file .sqlitedb in QGIS or ArcGIS? You can file convert? Database structure filesqlitedb is from DB Browser


texts: QGIS 2.4: copying of points to another point layer

texts:
in QGIS 2.4  it is not possible any more to copy points of a layer into another layer - I get the message:
" Fehler: Objektschreibfehler:
Objekterzeugungsfehler (OGR-Fehler: Attempt to write non-point (MULTIPOINT) geometry to point shapefile.) "

Is it a mistake of QGIS 2.4 ? In version 2.2 it was no problem. The two point-layers have the same structure (attributes)

My QGIS runs under ubuntu 14.04


texts: Accidentally drew polygon in WGS on top of UTM Map in QGIS?

texts: Is there any way to Change the polygon projections so it remains in place, but takes on the UTM coordinates?


texts:
With leaflet with reference to this sample: https://github.com/stuartmatthews/leaflet-geotiff
I was able to host a GeoTIFF file on the server and access it without tiling.
Can the same be achieved with Mapbox GL JS?
I have been looking into Mapbox GL JS documentation however could not find any information on how to access GeoTIFF hosted on a server.


texts: How to make a polyline using distance and direction?

texts:
QGIS newbie here. Given a table of values for distance (km) and direction (degrees), I want to create a polyline starting from a coordinate (lon, lat). I am going to use this to backtrack the origin of a particle in the ocean. My question is very similar to this:

How to draw a polygon given a distance and bearing in QGIS

The problem is that I do not know how to use PyQGIS.


texts: Export a list of invalid geometries using Python

texts:
Is there a way to export a list of layers containing invalid geometry (e.g. self-intersections) via Python?

Our GIS server contains over 40,000 layers, some of which have all sorts of errors that simply cannot be identified manually.

Is there some sort of OGR/GDAL command that can perform this?


texts: Esri leaflet layers turned off by default in map layer control box

texts:
I have a series of esri leaflet layers in a variable called "baseLayers" but want to initiate the map with all of them turned off by default inside the layer control box. I've tried playing with the L.control.layers method options but can't find anything that works. What I have currently is this:

L.control.layers(null, baseLayers, {collapsed:false}).addTo(map);

The 'addTo(map)' command is what I think causes all layers to show up as enabled in the layer control box by default, but if I remove that part of the code then the layers just don't show up at all in the control box. Any ideas?


texts: Creating a runtime visualization with CartoDb

texts:
I am trying to create a runtime visualization for a CartoDB map and I am quite the trouble with doing it. I have gone over the documentation and cannot get the code to work correctly at all unless it is a static vis. We are wanting to use the runtime code because our data is ever updating and we want the map to update with the code. Here is what I have thus far.

<link rel="stylesheet" href="http://libs.cartocdn.com/cartodb.js/v3/3.15/themes/css/cartodb.css" />

After writing this code and following along with the cartoDB documentation I don't get an output in the browser or any JS errors in the console. My hope is that I can receive support on how to correctly write the code that will create a runtime visualization of a map from our dataset.

I am not sure if what I am wanting is possible or not. We can create a static visualization however, when using the cartodbJS I am having no luck on creating a carto map from scratch using the above code.


texts: Adding feature class to Map in CURRENT project using ArcPy in ArcGIS Pro?

texts:
I want to add a feature class within a Geodatabase to a map in ArcGIS Pro (1.4.1) using arcpy. I am converting some python code from ArcGIS Desktop to ArcGIS Pro.

This issue is very similar to:

Adding newly created layer to Map in CURRENT project using ArcPy for ArcGIS Pro? and
Adding feature class to map (within specific group) using ArcPy with ArcGIS Pro?

However, instead of adding a feature class (or layer) to a map using the Python Window I am doing this from a  Script within a Toolbox. I have included a code snippet below.

lyrTest = r"C:\data\test.gdb\Layer1"
aprx = arcpy.mp.ArcGISProject("CURRENT")
aprxMap = aprx.listMaps("MainMap")[0]
lyrFile = arcpy.mp.LayerFile(lyrTest)
aprxMap.addLayer(lyrFile)

I get the following error message:

Traceback (most recent call last):
File "C:\data\AddLayers.py", line 5, in <module>
aprxMap.addLayer(lyrTest)
File "c:\users\user\appdata\local\programs\arcgis\pro\Resources\arcpy\arcpy\utils.py", line 191, in fn_
return fn(*args, **kw)
File "c:\users\user\appdata\local\programs\arcgis\pro\Resources\arcpy\arcpy\_mp.py", line 1048, in addLayer
return convertArcObjectToPythonObject(self._arc_object.addLayer(*gp_fixargs((add_layer_or_layerfile, add_position), True)))
ValueError: C:\data\test.gdb\Layer1
Failed to execute (AddLayers).


texts: Get Polygon or IGeometry from Envelope in ArcObjects

texts: I would like to create a Polygon Feature from the geometry of an IEnvelope in C# ArcObjects SDK 10. Is there an easy way how to do it?


texts:
I'm using QGIS and trying to clip a raster image with a valid polygon
(using GDAL), both in the same SRC (EPSG:32622). But the resulting image does not result in a valid raster. Does anyone know how to fix it?

The values of the input raster:

On the options, i have checked the option to keep the resolution of the input raster. (It's in portuguese, but it is only to give you the picture, everything else is the default value)

The result is an invalid raster:

I found it particularly interesting that the min value is positive and the max value is negative (according to what QGIS is telling me).


texts: To identify overlapping lines between 2 shape file using GDAL

texts:
I have 2 shape files in and I have to identify the overlapping lines between 2 shape lines. I have tried below code using overlap, It provided me non-overlapping lines as well.

while ((jsonFeature = layer.GetNextFeature()) != null)
        {
            var gdef = jsonFeature.GetGeometryRef();
            layer1.ResetReading();

            while ((jsonFeature1 = layer1.GetNextFeature()) != null)
            {

                gdef1 = jsonFeature.GetGeometryRef();

                isoverlapp = gdef1.Overlaps(gdef);

                if (isoverlapp == true)
                {
                    Console.Write("True");
                }

            }
   }


texts: When asked to select raster for elevation (in 3D configuration/QGIS) will be the DEM raster as downloaded or DEM raster re projected?

texts:
Do you need to re-project Aster GDEM data from USGS Explorer into a metric coordinate system in order to be able to view a map in 3d mode in QGIS?
 Or you just use the file input as it is downloaded from USGS?


texts:
I created a user, password and a role, granted access to the layer 'myfeature' only for that role. Now I want to test whether the access restriction works, but I cannot find any documentation about how the parameters must be set for user access in a ordinary WFS request.

however, neither:

http://mygeoserver/wfs?service=WFS&version=2.0.0&request=GetFeature&typeNames=myworkspace:myfeature&

nor

http://mygeoserver/wfs?service=WFS&version=2.0.0&request=GetFeature&typeNames=myworkspace:myfeature&user=myusernamer&password=mypassword&

work, unless I disable the access restriction.