texts:
I have a few MODIS products and as I was combining my imges I realised that some of the pixels around the coastline where abit off.Is that a standard issue and if so should I use reproject or should I georeference them on another image.And if I have to georeference those two images, should I use ArcGIS or ERDAS (I would like purely your opinion on the last question because I just wanted to know which one does this community prefer)

The Products that I am using are from MYD13Q1(AQUA) and MOD13Q1(TERRA) and the problem isn't just on the coastlines, it is also on every lake of the image.


texts:
Why is it useful to store data in unusual projections and datums?

I understand the value of projections as output constructs, because of distortion, etc: that much is justifiable.

However, I don't understand why, for instance, states use state plane projections in data. We have decimal accuracy: why not just store precise values in EPSG:4326? Is this entirely a remnant of pre-auto-reprojection days, or is there a value proposition I'm missing?

I'll restrict the scope of this question to vector data only to make it more concrete.


texts: ExportToJPEG - can't change resolution, always resets to default

texts:
In ArcGIS Pro (latest version), using the ExportToJPEG (or even ExportToPNG) arcpy function, I can't seem to export an image with my desired resolution. It always defaults to 96 dpi even when I specify 500. What am I doing wrong? Is this a bug?

This is the script I paste into the python console within ArcGIS Pro:

p = arcpy.mp.ArcGISProject("CURRENT")
print ("project loaded")
l = p.listLayouts("MyLayout")[0]
print ("layouts listed")
if not l.mapSeries is None:
 ms = l.mapSeries
 if ms.enabled:
   for pageNum in range(1, ms.pageCount + 1):
     ms.currentPageNumber = pageNum
     print("Exporting page {0} of {1}".format(str(ms.currentPageNumber), str(ms.pageCount)))
     l.exportToJPEG(r"C:\\Users\\Theo\\Dropbox\\site_plans\\Current_GR_site_plans\\GR2-GroundMovement-SP-_{0}".format(str(ms.currentPageNumber) + ".jpg", 500))

You can see at the end I specify 500 as the resolution. Even when I write- resolution = 500, it's the same. My images aren't coming out in my desired quality.

If I export a single page using the share ribbon, setting the dpi to 500, it works. But I'm automating this, so that's not a viable option for me.

Guides I've used:

https://community.esri.com/external-link.jspa?url=http%3A%2F%2Fpro.arcgis.com%2Fen%2Fpro-app%2Farcpy%2Fmapping%2Fmapframe-class.htm

https://community.esri.com/external-link.jspa?url=http%3A%2F%2Fpro.arcgis.com%2Fen%2Fpro-app%2Farcpy%2Fmapping%2Flayout-class.htm

similar issue to:

https://community.esri.com/thread/201455-exporttojpeg-resolution-problem

ExportToJPEG() arcpy script failure
(ArcMap)

neither helped me.


texts: Converting Euclidean Distance into meters using ArcGIS Pro

texts:
I used Euclidean Distance tool on a vector layer (school, point) and in the output raster layer, the distances are in weird unit and the distance ranges from 0.01 to 0.32.
How do I convert the distance into meters?
I need to find area within 250m, (and above 500m, 750, etc.) of the school.
I am using ArcGIS Pro


texts: I need to generate a unique id using field calculator in ArcMap. I am a student and I am having problems finding the script or getting auto increment to work or anything else.


texts:
I am trying to create a postgis query that would get all the points in an union of circles (possibly more than a thousand of circles with different radius). The circles would be created dynamically from a set of objects with a latitude, longitude and radius.

Is it possible to achieve this with postgis without perfomance issues?


texts:
The following code used to extract the data from a history file.
import osmium

class CounterHandler(osmium.SimpleHandler):
    def __init__(self):
        osmium.SimpleHandler.__init__(self)
        self.num_ways = 0
        self.elemtimeline = []

    def way(self, n):
        self.elemtimeline.append([n.id,
                                  n.version,
                                  n.visible,
                                  pd.Timestamp(n.timestamp),
                                  n.uid,
                                  n.changeset,
                                  n.deleted])
        self.num_ways += 1

h = CounterHandler()

h.apply_file("test.osm.pbf")
colnames = ['id', 'version', 'visible', 'ts', 'uid', 'chgset','deleted']

elements = pd.DataFrame(h.elemtimeline, columns=colnames)
elements.head(10)

The output table:

When this is checked on OSM https://www.openstreetmap.org/way/30046162 it shows that this is deleted over 9 years ago.

However in the table the deleted column shows False. Why is the conflict? Or is there a way to find the deleted objects from the history file?


texts:
I have a parcel layer and a table of property sales.  I have the two tables related (it is a one-to-many relationship) in a QGIS project so that when you click on a parcel you see all the related sales information.  Is there a way to query that related sales table and have it select the appropriate parcels directly in QGIS?

I have come up with a round-about way where I create a view in the Spatialite DB that the data is stored in but it is asking my end user a bit too much to be able to change the views parameters every time he needs a new sales map.


texts: Display custom fields/content in popup for existing ArcGIS Map API Service

texts:
I'm using the ArcGIS Map API "Historic_National_Boundaries_NEW (FeatureServer) > Year_2000" present over here :

API Link

It has two pre-defined fields 'OBJECTID' and 'NAME_2' (for each country) which are displayed onclick of a country in map. I'm using PopupTemplate() object to create the popup.

//on-click popup functionality starts
        var popupOptions = {
        markerSymbol: new SimpleMarkerSymbol("circle", 32, null,
          new Color([0, 0, 0, 0.25])),
        marginLeft: "20",
        marginTop: "20"
      };
       var popup = new Popup(popupOptions, domConstruct.create("div"));
        var popupTemplate = new PopupTemplate({
        title: "{NAME_2}",
        fieldInfos: [
          {
            fieldName: "NAME_2",
            visible: true,
            label: "Country: ",

          },
          {
            fieldName: "OBJECTID",
            visible: true,
            label: "ID: ",

          }

        ],
        showAttachments: true
      });

I need to know if these fields can be customized or new fields can be added. I checked through the documentation but was not able to find any info. It mentioned adding fields, media elements and plain text but not external content or any method to manipulate the fields.

If fields cannot be manipulated, how can I display custom content (possibly fetched from JSON/XML file) in the popup for each country?


texts: Clipping data to shapefile using Google Earth Engine?

texts:
I want to clip this fire data so it only displays results in my shapefile titled 'tapajos'. How do I do this?

var dataset = ee.ImageCollection('FIRMS').filter(
    ee.Filter.date('2000-01-01', '2018-08-10'));
var fires = dataset.select('T21');
var firesVis = {
  min: 325.0,
  max: 400.0,
  palette: ['red', 'orange', 'yellow'],
};
Map.setCenter(-119.086, 47.295, 6);
Map.addLayer(fires, firesVis, 'Fires');


texts:
I have a list of points showing the locations of some planned site investigations. Due to practical issues some of the locations were changed during the field work. The new locations was surveyed with a GPS. I've managed to merge the two files by point-name.

What I'm trying to do is to create a mark (point) showing the planned position and a line between the planned and actual position of the investigation to show the "offset".

Tried this with the geometry generator:

make_line(centroid($geometry),make_point( "183000 LETBANE_BOR_DKTM3_Easting" , "183000 LETBANE_BOR_DKTM3_Northing" ))

Seems like the "centroid.."-part is working, but the "make_point.." from coordinates is giving the wrong result:


texts:
I would like to add an on-map drawing function in my stand alone app, so I used QgsRubberBand to make it. I made a tool that emits points and then draw the QgsRubberBand. The problem is I don't find any solution to style it (applying a qml file would be perfect).
Here is the code I use to draw my QgsRubberBand:

    def update_largage_eau_canadaire(self, model):
        points = model.largage_eau_canadaire
        #points is an array of QgsPoint
        if len(points) > - 1:
            if self.largage_c_rubber_ is None:
                self.largage_c_rubber_ = QgsRubberBand(self.canvas_, False)

            self.largage_c_rubber_.setToGeometry(QgsGeometry. fromPolyline(points), None)
            self.largage_c_rubber_.setColor(QColor(255, 0, 0, 255))
            self.largage_c_rubber_.setWidth(2)

Here is the result:

And I would like do draw a styled rubberband that would look like this:


texts: How to draw downhill path with TileMill and KML?


texts: ModelBuilder Output Creates .lock file and Cannot Overwrite When Ran a Second Time

texts:
I am having trouble overwriting an output shapefile from ModelBuilder.

The model works great when ran the first time, but I get an '000210 : Cannot create output .' error if I re-run the model. The error states 'The output cannot be created. Potential reasons include data locking, an incorrect path, and limited access rights..'

I've noticed that there is a .lock file associated with the shapefile so I assume this is my problem. I've checked the overwrite outputs of geoprocessing options checkbox in the geoprocessing options.

My model is very simple, I am just using a multivalue iteration to buffer a line segment at different buffer widths and then clipping to a polygon. I am using the inline variable subsitution based off of the buffer value for my output shapefile name.

Anybody know how to fix this? or possible work arounds? I am using ArcMap 10.2


texts:
The goal is to check whether the user has added a point to a layer, and wait until then
Unfortunately, using QtGui.QApplication.processEvents() did not suffice.
time.sleep() inside a while loop on condition does not work, since QGIS will freeze its own Python engine with time.sleep
This is my best effort:
from qgis.core import QgsVectorLayer

# create layer
points_layer = QgsVectorLayer("Point", "temporary_points", "memory")
QgsProject.instance().addMapLayer(points_layer)

# get feature coint in points layer before editing
count_before_editing = points_layer.featureCount()
print("count_before_editing", count_before_editing)

# check wether the layer is editable
if not points_layer.isEditable():
    #make it editable
    points_layer.startEditing()
    # check if the requestet point feature has been set
    while points_layer.featureCount()<count_before_editing+1:
        print(points_layer.featureCount())
        # here i am missing the functionality to wait, while I manually add a point to the points layer
        #break is only here, to prevent your qgis from remaining in while
        break
    # stop editing and commit changes
    points_layer.commitChanges()


texts:
Ok, here is the situation. I am looking to get a custom arcobjects vb.net extension running correctly on ArcMap10.5. When loading nodes in an attachment viewer I receive the following error once Private Sub LoadNodes() is called.

System.TypeLoadException: Types extending from COM objects should override all methods of an interface implemented by the base COM class.

If I remove these lines from LoadNodes() I don't get the System.TypeLoadException: and the program runs through LoadNodes() but obviuously fails because the Attachment table is missing.

   ' Get relAttachments table and set selection table for it
    Dim relAttach As ESRI.ArcGIS.Carto.StandaloneTable = New cAttachStandAloneTable(mFWS.OpenTable(modDBNaming.GetQualifiedTableName("RelAttachments", mFWS)), "FCName")
    mRelAttachSel = relAttach

Somehow this issue is related to the newer environment because this issue does not occur in the old environment.

Current Environment: Visual Studio 13 .net framework 4.5 ArcMap and Arcobjects 10.5 (10.4 has the same issue)

Old Environment: Visual Studio 8 .net framework 3.5 ArcMap and ArcAobjects 10.2.2 and older

Here are the beginning steps of LoadNodes()

Private Sub LoadNodes()
    Try
        lblStatus.Text = "Populating the tree view with the database layers"
        Me.Cursor = Windows.Forms.Cursors.AppStarting

    ' Get relAttachments table and set selection table for it
    Dim relAttach As ESRI.ArcGIS.Carto.StandaloneTable = New cAttachStandAloneTable(mFWS.OpenTable(modDBNaming.GetQualifiedTableName("RelAttachments", mFWS)), "FCName")
    mRelAttachSel = relAttach

    ' Load relAttachments into dataset datatable
    PopulateRelAttachmentsTable(relAttach.Table)
    Dim relAttachDV As DataView = New DataView(mDataset.Tables("relAttachments"))

    ' Clear out reachnodes hashtable
    mReachNodes.Clear()

    ' Run through any selected features and store in table
    If mStandAlone = False Then
        PopulateSelectionTable()
    End If
    Dim selDV As DataView = New DataView(mDataset.Tables("Selection"))

    ' Clear the tree
    tvStream.Nodes.Clear()

Code from frmAttachViewer

' Load the dictionaries
            LoadDictionaries()
            MessageBox.Show("Loaded Dictionaries")

            ' Load the featureclass names into cboFeatureclass
            For Each de As DictionaryEntry In mLayerTab
                cboFeatureClass.Items.Add(de.Key)
            Next

            ' Load the tree view nodes
            LoadNodes()

            ' "Open" the reaches that have been selected (or features in them were selected)
            mSuppressSelChanges = False
            If mStandAlone = False Then
                MapSelectionChanged()
            End If

        Catch ex As Exception
            StreamAnalyst_AGS10.ErrorHandler.ShowErrors(ex, "Warn", "An error occurred while loading the form")
            log.Warn("An error occurred while loading the form", ex)


texts:
I'm trying to convert Sentinel2 jp2 image to tiff using gdal_translate. When running the code:

Miguel$ gdalinfo /Users/Miguel/Desktop/S2_b1.jp2
Driver: JP2ECW/ERDAS JPEG2000 (SDK 5.2)
Files: /Users/Miguel/Desktop/S2_b1.jp2
Size is 1830, 1830
Coordinate System is:
PROJCS["WGS 84 / UTM zone 33N",
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AXIS["Latitude",NORTH],
    AXIS["Longitude",EAST],
    AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",15],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
    AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","32633"]]
Origin = (600000.000000000000000,7000020.000000000000000)
Pixel Size = (60.000000000000000,-60.000000000000000)
Metadata:
ALL_COMMENTS=Kakadu-v7.4
Kdu-Layer-Info: log_2{Delta-D(squared-error)/Delta-L(bytes)}, L(bytes)
-192.0,  3.8e+06

COLORSPACE=GREYSCALE
COMPRESSION_RATE_TARGET=3
Corner Coordinates:
Upper Left  (  600000.000, 7000020.000) ( 16d58'56.33"E, 63d 6'56.41"N)
Lower Left  (  600000.000, 6890220.000) ( 16d55' 3.52"E, 62d 7'50.52"N)
Upper Right (  709800.000, 7000020.000) ( 19d 9'13.82"E, 63d 4' 7.06"N)
Lower Right (  709800.000, 6890220.000) ( 19d 1' 7.24"E, 62d 5' 8.11"N)
Center      (  654900.000, 6945120.000) ( 18d 1' 5.23"E, 62d36'15.28"N)
Band 1 Block=256x256 Type=UInt16, ColorInterp=Gray
Description = Grayscale
Overviews: 915x915, 457x457, 228x228
Image Structure Metadata:
NBITS=15

ns27:~ Miguel$ gdal_translate /Users/Miguel/Desktop/S2_b1.jp2       /Users/Miguel/Desktop/S2_b1tif.tif
Input file size is 1830, 1830
0

The file is created but can´t be opened.

Any suggestion about how to proceed?


texts: Printing the values on the trendline in the time series chart to the console at GEE

texts:
How can I print the values above the trendline in the time series chart to the console? For example, I want to print the first and last number on the trendline to the console.
Code Link:
https://code.earthengine.google.com/?scriptPath=users%2Formserdar%2Fdefault%3Agis_stack
var startTime = '2000-01-01';
var endTime = '2016-01-01';

// x The polygon that describes the plot.
var plot = /* color: #d63000 */ee.Geometry.Polygon(
        [[[28.51818846650565, 38.69338161384325],
          [28.518188466504512, 38.6927510390996],
          [28.5189930726267, 38.69275103633082],
          [28.518993072625108, 38.693381611074464],
          [28.51818846650565, 38.69338161384325]]]);

//Add Map layers

var dot = ui.Map.Layer(plot, {color: 'FF0000'});
  Map.layers().set(1, dot);
  Map.centerObject(plot);
  Map.setOptions('satellite');

//-----------------------------CHARTS-------------------------------

//Landsat 7 32-Day NDVI Composite
var landsat7MonthlyNdvi = ee.ImageCollection('LANDSAT/LE7_L1T_32DAY_NDVI').filterDate(startTime, endTime).select('NDVI');
var landsat7MonthlyTimeSeries = Chart.image.series(landsat7MonthlyNdvi, plot, ee.Reducer.mean(), 30);
landsat7MonthlyTimeSeries = landsat7MonthlyTimeSeries
    .setOptions({
        title: 'Landsat 7 32-Day Monthly NDVI Composite',
        hAxis: {title: 'Date', gridlines: {count: 10}},
        vAxis: {title: 'NDVI',viewWindowMode: 'explicit', viewWindow: {max: 1,min: -0.25,},gridlines: {count: 5,}},
        interpolateNulls: true,
        lineWidth: 1,
    pointSize: 1,
    crosshair: { trigger: "both", orientation: "both", color: 'green' },
    trendlines: { 0: {type:'linear', showR2: true, lineWidth:2, pointSize: 0, color:'red', visibleInLegend: true}}
    });

print(landsat7MonthlyTimeSeries, landsat7MonthlyNdvi);


texts: Is it possible to display two categorical symbolizes for two different attributes in one layer alongside each other in ArcMAP 10.1?

texts:
I have a road layer that I would like to display attribute_1 = "speed" on top of the line segment and attribute_2 = "traffic_volume" below the line segment.
For each attribute I would like to setup speed bins 0-30, 30-40, 50-60, 60+ and for volume bins: 0-200, 200-1000, 1000-3000, 5000-10000, 10000+
As a rough Schematic, red represents speed, and black represents volumes and the road should run between those two symbols:

Question: I know how to display each attribute as symbol; however, I'm not sure what's the best way to display both symbologies so that they don't overlap
Note that I'm open to using ArcPY


texts: Mapnik XML mapfile projection of multiply layers

texts:
How to apply a multiply layer with Mapnik from EPSG:4326 TIFFs to EPSG:3875 projected tiles?

I have two TIFF files, one containing the color-relief and the other the slope-shade (generated with GDAL, which is perfectly demonstrated here). The slope-shade is going to be blended with the color relief image, using the multiply blend mode. Both TIFFs are in EPSG:4326.

With Mapnik, to include the input projection (EPSG:4326) and the output projection (EPSG:3875) You have to add the srs to the Map and both Layer elements. This is done in the example file below. Here is an explanation by Frederik Ramm about the approach.

The color-relief map is projected correctly, but the slope-shade is not multiplied. How can this be explained and how should the projection be modified to actually multiply the slope-shade layer with the color-relief layer?

<Map srs="+init=epsg:3857">

  <Style name="color relief style" >
    <Rule>
      <RasterSymbolizer mode="normal" />
    </Rule>
  </Style>
  <Style name="slopeshade style">
    <Rule>
      <RasterSymbolizer opacity="0.4" mode="multiply" />
    </Rule>
  </Style>

  <Layer name="color relief" srs="+init=epsg:4326">
    <StyleName>color relief style</StyleName>
    <Datasource>
      <Parameter name="type">gdal</Parameter>
      <Parameter name="file">/path_to_dir/map_color.tif</Parameter>
    </Datasource>
  </Layer>
  <Layer name="slopeshade" srs="+init=epsg:4326">
    <StyleName>slopeshade style</StyleName>
    <Datasource>
      <Parameter name="type">gdal</Parameter>
      <Parameter name="file">/path_to_dir/map_slopeshade.tif</Parameter>
    </Datasource>
  </Layer>
</Map>

Note: When I transform the TIFFs to EPSG:3857 and replace the Layers srs with srs="+init=epsg:3857" in the mapfile, the slopeshade is multiplied correctly. But, I don't want to transform my TIFFs for other reasons.


texts: How to Create Query Status Table Postgres

texts:
I am currently a user of a Postgres database (I am not the DBA). Within my sandbox I run several business critical scripts every day and currently track their status by first in a log text file on my desktop while running the scripts in a batch file. I am trying to create a table in my sandbox that can be used for reporting and analysis that includes the status of the scripts as well as if any of them have failed to run properly or not.
Such as:

----------------------------------------+-----------
 Beginning Script-A                     | Timestamp
 Script-A Completed Successfully        | Timestamp
 Beginning Script-B                     | Timestamp
 Script-B failed to Run Successfully    | Timestamp
----------------------------------------+-----------

How could I do this?


texts: How to exclude water during contour generation?

texts:
I recently acquired a LiDAR data set that included a 1 meter DEM. I generated 1, 2 and 5m contour layers but I am seeing a lot of contours within water bodies. I suspect this is either caused by wave movement or in shallow areas where rocks/debris are present. I generated the contours using ArcGIS Desktop 10.3.1 3D Analyst.

What are my options for removing these contours? Should this elevation issue be addressed at the Lidar/DEM level? I have a shapefile that delimits the land from water that I could use for a break line.


texts: Using Geoserver and PostgreSQL to build a web mapping tool?

texts:
I've been an ESRI user and a Carto user and like both tools but I am now hoping to build my own web-mapping platform. The issue is, I have no web-development experience other than the little I've taught myself. I want my mapping platform to be able to run SQL queries and to store data safely (some datasets I only want to make available via a link, and not make the data shareable/scrapeable).

Am I even looking in the right direction?

Can this GeoServer + PostgreSQL mix get me what I want?


texts: How can i plot this piece of code with ggplot2?

texts:
I have these vectors of the lakes in New Guinea and the coastline and  i want to plot it using ggplot2 exactly as it follows but with a more cool appearance.

co<-coastline.shp
nglakes<-lakes.shp
plot(co,col="yellow")
plot(nglakes,col="blue",add=T)
legend("topleft",legend=c("LAND","LAKE"),title="LineTypes",col=c("yellow","blue"),lty=1,cex=0.8)


texts: OpenLayers 3 custom control doesn't get the map object

texts:
I've been playing around with extending controls in OpenLayers 3. I've built a control which grabs a legend image from a mapserver with a GetLegendGraphics request. However, when I add the control to the map, the map object doesn't get passed to the control. I think the inheritance works just fine, because the control does have a getMap() function, but it returns undefined. Does anyone have an idea what am I doing wrong?

Code:

var ol3_legend = function(e) {
    var options = e || {};
    var wmsVersion = options.wmsVersion || '1.3.0';
    var format = options.format || 'image/png';
    var legendP = document.createElement('p');
    legendP.innerHTML = 'Legend:';
    var legendDiv = document.createElement('div');
    legendDiv.className = options.class + ' ol-unselectable';
    legendDiv.appendChild(legendP);
    var layers = this.getMap().getLayers().getArray();
    for (var i=0;i<layers.length;i++) {
        if (layers[i].get('showLegend') === true) {
            try {
                var url = layers[i].getSource().getUrls()[0];
            }
            catch(err) {
                var url = layers[i].getSource().getUrl();
            }
            var legendImg = document.createElement('img');
            legendImg.src = url + '&version=' + wmsVersion + '&service=WMS&request=GetLegendGraphic&sld_version=1.1.0&layer=' + layers[i].getSource().getParams().layers + '&format=' + format;
            legendDiv.appendChild(legendImg);
        }
    }
    ol.control.Control.call(this, {
        element: legendDiv
    });
}
ol.inherits(ol3_legend, ol.control.Control);

ol3_map.addControl(new ol3_legend({
    class: 'ol_legend'
}));

Output:

Uncaught TypeError: Cannot read property 'getLayers' of undefined

Workaround (just the modified lines):

var layers = options.map.getLayers().getArray();

ol3_map.addControl(new ol3_legend({
    map: ol3_map,
    class: 'ol_legend'
}));

The problem is, that I have to provide the map object, which makes restrictions if I would like to reuse the code later.

Working fiddle.


texts: Converting from .las to .ply and vice versa in Python?

texts:
I am trying to convert a .las/.laz file into .ply file, read the .ply file, do some processing on it and then convert it back to .las/.laz file.

Is there any Python implementation for this?

For reading a .las file, I use the below code:

inFile = laspy.file.File(folder + '.las', mode='r')

I have tried searching for next steps and found this but unable to use it as I want.


texts:
I am new to GIS and I want to compare two DEM raster .tiff files, i.e. I want to compute differences between height values etc.

How can I determine if the two tiff files are comparable?

I have read that SRTM is measured relative to sea level, here. So is it measured relative "mean sea level"? Or is that equivalent?

The other tiles were given originally based on the height reference 'DE_DHHN2016_NH' with crs 'ETRS89_UTM32' and geoid 'DE_AdV_GCG2016_QGH', documentation is here.

The SRTM tile dgalinfo gives me:

Driver: SRTMHGT/SRTMHGT File Format
Files: N50E007.hgt
Size is 3601, 3601
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]
Origin = (6.999861111111112,51.000138888888891)
Pixel Size = (0.000277777777778,-0.000277777777778)
Metadata:
  AREA_OR_POINT=Point
Corner Coordinates:
Upper Left  (   6.9998611,  51.0001389) (  6d59'59.50"E, 51d 0' 0.50"N)
Lower Left  (   6.9998611,  49.9998611) (  6d59'59.50"E, 49d59'59.50"N)
Upper Right (   8.0001389,  51.0001389) (  8d 0' 0.50"E, 51d 0' 0.50"N)
Lower Right (   8.0001389,  49.9998611) (  8d 0' 0.50"E, 49d59'59.50"N)
Center      (   7.5000000,  50.5000000) (  7d30' 0.00"E, 50d30' 0.00"N)
Band 1 Block=3601x1 Type=Int16, ColorInterp=Undefined
  NoData Value=-32768
  Unit Type: m

I merged the other tiles to get one large raster. Here are the details given by gdalinfo:

Driver: GTiff/GeoTIFF
Files: bonn_dsmEPSG_4647.tif
Size is 14000, 18000
Coordinate System is:
PROJCS["ETRS89 / UTM zone 32N (zE-N)",
    GEOGCS["ETRS89",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6258"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4258"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",32500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","4647"]]
Origin = (32359999.500000000000000,5627999.500000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (32359999.500, 5627999.500) (  7d 0'49.80"E, 50d47'12.56"N)
Lower Left  (32359999.500, 5609999.500) (  7d 1'14.38"E, 50d37'30.16"N)
Upper Right (32373999.500, 5627999.500) (  7d12'44.43"E, 50d47'24.12"N)
Lower Right (32373999.500, 5609999.500) (  7d13' 6.55"E, 50d37'41.66"N)
Center      (32366999.500, 5618999.500) (  7d 6'58.80"E, 50d42'27.28"N)
Band 1 Block=14000x1 Type=Float32, ColorInterp=Gray


texts: Interpolation of contours (from one "main-contour-shapefile")

texts:
I calculated, that the best and smallest possible äquidistance (equidistant) for my map is 30m. After I generalized my DHM, I extracted contours of a interval of 30m and smoothed them. Now I have my "main-contours".

But I also want contours with 10m äquidistance (equidistant), and I need to interpolate it from my main-contours (because my DHM is not accurate enough to extract 10m-contours directly).

EDIT: The 10m-contours are just for optical reasons (counting) and don't have to be accurate. The accuracy of the 30 m-contours is fine enough!

Is that possible with Grass, Saga or Qgis?

EDIT: The 30m-contours are accurate, (vertical error less then 15m), the 10m-contours have errors so I don't want them! :)


texts: How to resolve error when launching OSeo4W Shell: Windows cannot find .bat file?

texts:
I downloaded and installed OSGeo4W to use ogr2ogr, but the the program OSGeo4W will not run. On launching the shell, I receive an error stating the .bat file cannot be found in specified directory (image included) even though the .bat file is in the specified directory.

When I run OSGeo4W as an administrator the shell will open with the error of filename, directory name, or volume label syntax is incorrect (image included).

I have uninstalled and reinstalled OSGeo4W in an attempt to resolve the error; however, I don't think it is an installation issue since I was using it yesterday without producing the described error.

The only thing that I've done since yesterday that could have effected OSGeo4W is install Cygwin on my computer, which runs Windows, to use batch.


texts: Compare spatial cluster/agglomeration patterns for polygon data

texts:
I have agricultural data in form of polygons which I would like to test for spatial clusterings/spatial agglomerations.

All in all I have about 40 variables which I can aggregate and standardize in different ways. One way of standardization could be for example to calculate production values per capita within each polygon. Another way could be to calculate production values per ha within each polygon.

All ways of standardization and aggregation produce different maps with different spatial patterns: clusterings and no-clusterings.
So as a base for my later analysis I wan't to identify such aggregation/standardization combinations that produce  strong spatial clusterings. Therefore I would need to compare the different results from aggregation and standardization.

Of course I could do this looking manually at the maps (see example below). But this is quite subjective and only in some case you can make clear distinctions. Imagine doing this for 40 variables and let’s say 8 possible ways of data preparation… So I would prefer to use some objective measurement i.e. spatial statistics.

I use R and Arc GIS. Has anyone an idea how to implement such an analysis?

The examples below show Banana Production once without standardization and once standardized per capita. they look very similar, but which one is more spatially clustered?


texts: Calculate geometry on some features have same X and Y for start and end?

texts:
I did "Calculate Geometry" on a line shapefile to obtain the coordinates X and Y for the start and the end of all the lines of the shape. But some lines resulted with the same X on the start and on the end, and the same Y on the start and on the end.

If the line is crossed by other line, the X and Y that results is the X and Y of the point where the line is crossed. So X,Y are the same because it's a point.
My question now is why Calculate Geometry does that?

This is one of the affected lines:

Thsi is another line and its properties


texts: Curved lines turn to straight lines after sDNA preparation

texts:
I am an architect, sometimes I need to use sDNA to analyze networks mapped in AutoCAD. But I met a problem with curved lines, i.e. circle, or polylines.

After importing .dwg file into .gdb, and using 'prepare' function of sDNA toolbox, all of the curved lines turns into straight lines with only starting and ending points. Is there any way to solve this problem?


texts: Changing projection of Google Earth Engine code editor display

texts:
I would like to change the projection on the code editor display window from EPSG:3857 (maps Mercator) to something else like Lambert conformal conic.

If this is not possible on code editor display window, what are my options to display Earth Engine images on a different projection?

I tried to look for tips on the projections documentation page but it can't find how to change the projection of the Map object in code editor.


texts: Drawing box arround text in OpenLayers 3?

texts:
I need text labels with surrounding a box and a specific fill color. Actually I have the following function which only shows the text but I don't know how to draw a box around the text. Text length is always different, so it must consider the text length. Can anybody show me how to implement it?

function showLabel(mapPos) {
  var textStop = replaceAll(mapPos['text'],"|","");
  //textStop = replaceAll(textStop," ","&nbsp;");
  //textStop = replaceAll(textStop,"-","&ndash;");
  var lblFeature = new ol.Feature({
    geometry: new ol.geom.Point(ol.proj.transform([mapPos['lng'],mapPos['lat']],'EPSG:4326','EPSG:3857'))
  });
  var textPos = new ol.style.Text({
    textAlign: 'left',
    textBaseline: 'middle',
    //font: 'bold 16px arial',
    font: '16px arial',
    text: textStop,
    fill: new ol.style.Fill({color: '#000000'}),
    //stroke has bad performance
    //stroke: new ol.style.Stroke({color: mapPos['color'], width: 1}),
    offsetX: 16,
    offsetY: -16,
    rotation: 0
  });
  //TODO: Draw a box arround the text
  var lblStyle = new ol.style.Style({
    text: textPos
  });
  lblFeature.setStyle(lblStyle);
  return lblFeature;
} //End showLabel()

So far I do not get any error in the console.log. I guess it must be done near my TODO Label. mapPos is an array where in mapPos['text'] is text available.

In version 2 I was able to do it like this function:

var markerFlag = new Array();
function showPosition(mapPos) {
  var size = new OpenLayers.Size(24,24);
  var icon = new OpenLayers.Icon('../images/map24/'+mapPos['symbolid']+'.gif',size);
  var myLonLat = new OpenLayers.LonLat(mapPos['lng'],mapPos['lat']).transform(
    proj4326, //Transform from WGS 1984
    projmerc //To Spherical Mercator Projection
  );
  var textPos = replaceAll(mapPos['text'],"|","");
  textPos = replaceAll(textPos," ","&nbsp;");
  textPos = replaceAll(textPos,"-","&ndash;");
  var myOptionsText = {
    mouseOver: false,
    labelOffsetTop: '-58px',
    labelOffsetLeft: '22px',
    backGroundColor: mapPos['color']
  };
  layer_position.addMarker(new OpenLayers.Marker(myLonLat,icon.clone(),myOptionsText));
  if (mapPos['text']) {
    layer_banner.addMarker(new OpenLayers.Marker.Label(myLonLat,icon.clone(),textPos,myOptionsText));
  }
} //End showPosition()


texts:
I'm newbie on this. I got some GIS data on Arcview but I want to develop a Web based GIS, I've found that OpenLayers is pretty easy to me because I'm Javascript/jQuery developer, the issue is that I don't have a clue how to export the data from Arcview to get it work with OpenLayers.

It would be nice if I can change some data over OpenLayers and still maintain those updates on the original data on Arcview.


texts:
What i have:
I've installed Qgis 2.99.0 (dev)

How do I:
Install base map for this version. The openlayers plugin isn't available.


texts:
I have run into a problem that I see has been posted previously, but none of the posted solutions seem to solve my problem.

I have:

A shapefile (with geometry) with a single row for many counties.
A CSV file (without geometry) that has multiple rows (representative of different time points) for each county.

Each of these files has an identifier in common “CQID”, which is unique to the county. As such, the shapefile has 1 of each CQID. The CSV has nine of each of the CQID’s (since I have nine time points).

I would like to join the CSV file into the shapefile. So far, I have tried and failed, to solve this problem doing the following:

Using the Join Feature in the Layer Property of the shapefile.
Creating a “Relation” of the two files, via Project properties

I expect to look at the attribute table of my shapefile, and see that where I had a single line for each county, there are now nine lines for each.

How do I achieve this is QGIS?

Alternatively, is it easier to do so in ArcMap? Any good tutorials for that?


texts: Set attribute length when exporting features with st_write

texts:
I have sf features that I have created/manipulated in R that I need to export to shapefile (my client only accepts shapefile format so I cannot use .gpkg or other formats). I have used st_write to successfully save the features as a new file, however, the length of each attribute is the maximum length possible. How can I set the length for each attribute? The help file indicates that this is probably possible but I can't quite figure out how. I'm guessing that is it is in dataset_options = NULL, layer_options = NULL.

Example using sample dataset

nc <- st_read(system.file("shape/nc.shp", package="sf"))
st_write(nc[,c("AREA", "NAME", "CNTY_ID")], "nc_test2.shp",name=c("area", "name", "ID"))

All the numbers are considered as float/double values with the maximum length and precision and the text values have a length of 80. How can I set the length when doing the export? In SQL, the statement would be something like this:

ALTER COLUMN (NAME TYPE CHAR(20),
        AREA numeric(5, 3),
        CNTY_ID int(5))


texts:
I'm using Winpython 3.6.7 and I want to create a standalone application using QGIS 3.4.11. When I append the QGIS path using the sys module, nothing happens.

import sys
sys.path.append("C:\OSGeo4W64\apps\Python37\lib")
import qgis.core
>>> ModuleNotFoundError: No module named 'qgis

Previously when I wanted to work with PyQGIS 2.x (built-in Python), I created a batch file to familiarize my IDE with QGIS and set PATH and PYTHONPATH variables to QGIS \bin and \python. e.g. Something like this for Pycharm IDE:

SET OSGEO4W_ROOT=C:\OSGeo4W64
SET QGISNAME=qgis-ltr
SET QGIS=%OSGEO4W_ROOT%\apps\%QGISNAME%
SET QGIS_PREFIX_PATH=%QGIS%
SET PYCHARM="C:\Program Files\JetBrains\PyCharm Community Edition 2017.1.4\bin\pycharm64.exe"

CALL %OSGEO4W_ROOT%\bin\o4w_env.bat
CALL %OSGEO4W_ROOT%\bin\gdall6.bat

SET PATH=%PATH%;%QGIS%\bin
SET PYTHONPATH=%QGIS%\python;%PYTHONPATH%

start "PyCharm aware of QGIS" /B %PYCHARM% %*

But now that I have to use Winpython and I want to import QGIS modules to it, this doesn't work. Basically, how we can integrate different python distributions to use QGIS library while enjoying full functionality? (To add some QGIS functionalities to a previously developed standalone application, for instance)

And yes, I know it's counter-intuitive to use another Python distribution (WinPython 3.6.7) other than the Python which is shipped with QGIS. But in this case, I have to. I would be delighted to see your suggestions.


texts: Rasterize a geometry using gdal.RasterizeLayer()

texts:
I am writing a processing script using QGIS 3.4. I need to rasterize a polygon using gdal.RasterizeLayer. The problem is the polygon is a list of QgsGeometry objects and not a layer.

n = (2*radiusInt)+1
epsg = bldgs_lyr.sourceCrs().authid()[5:]
srs = osr.SpatialReference()
srs.ImportFromEPSG(int(epsg))
geotransform = (sourceX-radius, 1, 0, sourceY+radius, 0, -1)

output_rast = gdal.GetDriverByName('GTiff').Create(raster_uri, n, n, 1, gdal.GDT_Float32)
output_rast.SetGeoTransform(geotransform)
output_rast.SetProjection(srs.ExportToWkt())

mem_layer = QgsVectorLayer("Polygon?crs="+bldgs_lyr.sourceCrs().authid(), 'polygon', 'memory')
prov = mem_layer.dataProvider()
for geom in geometries:
    feat = QgsFeature()
    feat.setGeometry(geom)
    prov.addFeatures([feat])

gdal.RasterizeLayer(output_rast, [1], mem_layer, burn_values=1)

As you can see, I've tried creating a vector layer to use as input. But if I run this, QGIS will just crash.


texts: Differentiating Cloud and Snow in unsupervised classification of mountaneous terrain?

texts:
I am doing unsupervised classification of a mountaneous terrain and I have used 500 classes to classify the land use.

In my result, cloud and the snow are categorized in the same class.

How do I differentiate them?

I want to eliminate the cloud but not the snow.


texts: Spectrum Spatial Analyst v12 issues on upgrade

texts:
After I upgraded spectrum ver 11 to 12 and spatial analyst 8 to 12 - my maps open but I don't get the layer control? In the SSA admin console I get a box to choose a color which doesn't go away and all the configs give errors.

Has anyone else tried and seen the same issues -I want to get this running during the weekend and so I can't hitup Pitney Bowes support.

Spatial Manager

The map view should look like (this is from our user guide)


texts: IDE Problem ArcGIS Runtime SDK for Android

texts:
I'm trying to develop an app with ArcGIS Runtime SDK for Android. Honestly, I'm not much experienced.

Just looking for the best IDE choise.

Eclipse works just fine with the Runtime SDK. But has no library documentation in the IDE.

IntelliJ IDEA 12+ is a good option for Android Development. But I couldn't use the Runtime SDK (v10.1.1) in it.

I want to import the documentation to the one of these IDEs. Such as :


texts: is possible creating Geowebcache tiles offline

texts:
When I start geowebcache seed mechanism, my server CPU is working %120 persentage. And I can not use WMS service on my openlayers application because it is working very slow. Server sources are consuming by geowebcache.

Can I create tiles offline via QGIS or another application.

Geowebcache cache directory is like this:

cache/topp_states/EPSG_2163_00/0_0/00_00.png
cache/topp_states/EPSG_2163_00/0_0/00_01.png
cache/topp_states/EPSG_2163_00/0_0/01_00.png
cache/topp_states/EPSG_2163_00/0_0/01_01.png


texts:
How do I multiply an image collection by a number?
I would like to multiply the image collection by 10 on Line 41.  I have named the variable and stated the variable to be multiplied on that line.  I have also adjusted the code to work with the variable I have defined.
The study area is currently an arbitrary small polygon drawn in GEE.  You can ignore the other study area assets.  The expected output will be in signed int16 form.
https://code.earthengine.google.com/e25ab68f0d7756032c7782c6768b56db


texts: Switching road names from uppercase to mixed case in ArcMap?

texts:
Our county road data has road names for all roads in the county. In the attribute name field, there are some entries with all uppercase for the road name and some with mixed case. I would like to go about switching the road name attributes all to mixed case just for uniformity throughout the county data.

How can this be done?

I know I could assign someone to manually edit every record into Mixed Case, but I'd like to try to perform this as one single efficient process if possible.

I've looked into Field Calculator and it looks like there could be an answer there, but I'm having trouble identifying which of the processes would work. I don't mind exporting it and working the table in Excel and then rejoining it after it's cleaned if that's truly the best way. We have 33,906 road segments.

I've attached some screenshots as examples.


texts: How to create a Legend with Points. Same size in Legend and Map. Leaflet

texts:
I am trying to create a Legend in Leaflet that contains the circles in the map.
But, the size in the legend should the same size that it is in the map.
My code is next.
function getRadius(pop) {
  var maxSymbolsize = 20; // maximum symbol size
  var maxValue = 37393129; // highest population value in the dataset
  r = maxSymbolsize * Math.sqrt(pop/maxValue); // proportional by area
  return r;
}

function getColor(d) {
return d > 10000000  ? '#d7301f' :
       d > 5000000   ? '#fc8d59' :
       d > 1000000   ? '#fdcc8a' :
                       '#fef0d9' ;
}

var legendcolor = L.control({position: 'bottomright'});
 legendcolor.onAdd = function (map) {
    var div = L.DomUtil.create('div', 'info legend'),
        grades = [0, 1000000, 5000000, 10000000],
        labels = [];

    for (var i = 0; i < grades.length; i++) {
        div.innerHTML +=
        ' ' +
        grades[i] + (grades[i + 1] ? '&ndash;' + grades[i + 1] + '
' : '+');
          }    return div;
   };
  legendcolor.addTo(map);

So, it gives me as a result:

As you see, it has circles but I want the sizes both the same in legend and map.
Here you see the CSS if can be involved in the solution. Where can I include the function GetRadius to make possible in legend.
.legend .circle {
  border-radius: 50%;
  width: 15px;
  height: 15px;
  margin-top: 8px;
}


texts: ApplySymbologyFromLayer not including all values in my domain

texts:
I am using a combination of ArcMap and ArcPy to stage data for sharing as a hosted Feature Service on AGOL.

The data will be used in a data-driven pages type of way, where each feature will be reviewed and have a 'ReviewStatus' attribute updated before moving to the next feature for the same review.

The way I have it set up, every feature has a default value in my tracking field of 'Awaiting Review'. As the feature is reviewed, a field is updated to indicate whether the changes were 'Accepted', 'Rejected', or 'Marked for Revisit'.

To facilitate this, I have created and applied a domain consisting of only those four values discussed above and created a layer file symbolizing each value as a different color.

When I have all my layers loaded on my desktop session, I am attempting to loop through each layer and apply the symbology from my staged layer file. The issue is...each feature has a default value of 'Awaiting Review', so before any review has happened, there are no features with any other values in my domain. When I use the ApplySymbologyFromLayer tool, it only brings the symbology for 'Awaiting Review'. It ignored 'Accepted', 'Rejected', and 'Marked for Revisit'.

I have tried messing with FeatureTemplates to no avail.

The manual workflow that HAS been successful is to simply import the layer file as symbology through the layer properties. When I close the properties dialog, everything works as expected. It is trying this pythonically that fails.

In my routine, I have also tried the layer method addAllValues() AFTER applying symbology and before. I even tried to reload the symbology after adding all values.

I know arcpy.mapping can be a bit clunky sometimes and often needs some inelegant workaround to accomplish a task, but am I missing something? Is there a way to include the symbology for ALL values in a layer file even if the dataset doesn't yet contain features with all those values?


texts: what dictates the drawing order of overlapping features with in the same layer/feature class in ArcGIS?

texts:
I know that in ArcGIS drawing order is controlled by layers, and if I want to make some features be drawn on top of other features I can use definition queries to separate them into different layers and position one on top of the other.

But I was wondering (mostly out of curiosity) what dictates the drawing order of the 1000 overlapping polygons in my feature class, and if there's a way that I can control what features are displayed on top within one layer?


texts: Calculating the area of the polygon segments within multiple overlapping buffers

texts:
Let me try to explain my problem:

I have a layer of polygons
A layer of overlapping buffers
I want to calculate the area per buffer of the containing polygons.

The problem is that the buffers are overlapping so the clip function does not understand that one polygon can belong to multiple buffers and at the same type each buffer should only have the area of the polygon within itself.
So far, I have tried to clip the polygons within the buffers, and then summarize within the areas, but somehow it does not get the correct answer.
I have also tried the weighted summarize within, but is not able to finish successfully (error writing output file)
I am thinking the only way is in notebook to write a for loop that takes buffer per buffer and summarizes the part of the polygon within. But I am not sure how to write the code for this.
Any pointers?


texts: Are there downloadable samples of Geoprocessing tools that can be published as Geoprocessing Service tasks?

texts: I am looking for examples of Geoprocessing Services which contain a tool/script for download.


texts:
I am trying to display an image in QGIS cloud when hovering over a point. I went through this tutorial:
https://blog.sourcepole.ch/2020/10/19/image-tooltip-with-qgis-cloud/
But I get this error:

This is the py. file that the site provides:
# Import all nesseccary Modules and Classes from PyQt5
from PyQt5.QtWidgets import QTextEdit,  QToolButton,  QFileDialog,  QLineEdit,  QLabel
import base64

# This is the main function which will be connected with the QGIS layer
def form_open( dialog, layer, featureid ):

# Define some global variables
    global name_field
    global image_field
    global image_preview
    global tool_button
    global file_name

    name_field = dialog.findChild( QLineEdit, "name" )
    file_name = dialog.findChild( QLineEdit, "filename" )
    image_field = dialog.findChild( QTextEdit, "image" )
    image_preview = dialog.findChild( QTextEdit, "imagePreview" )

# Find the QToolButton Object to connect the conversion function.
    tool_button = dialog.findChild( QToolButton, "tool_button" )

# When the ToolButton is clicked run the method get_file_name()
    tool_button.clicked.connect(get_file_name)

# set the imageField Widget invisible for a nicer GUI. It is not nesseccary to see the data
    image_field.setVisible(False)

#Make the image preview visible with the HTML image tag

    image_preview.setText('' % (image_field.toPlainText()))

# Open FileDialog and convert the binary image to base64
def get_file_name():
    file_path,  filter = QFileDialog.getOpenFileName(None,  'Open Image File',  '',  '(Images *.jpg *.JPG)')

    if len(file_path) > 0:
        file_name.setText(file_path)
        encoded_image = ''

# Read the image file and proceed the encoding
        encoded_image = base64.b64encode(open(file_path, 'rb').read()).decode('utf-8')

# Set the appropriate fields for storing in the QGIS layer
        image_field.setPlainText(encoded_image)

# Set HTML img tag around the base64 image
        image_preview.setText('' % (str(encoded_image)))

This is the result:

Not sure where I am going wrong here.


texts: Getting all populated places per state when number of places over 6000

texts:
I am able to query for the number of places.  The problem becomes when there are over 6000 places, the system gives me an error message that the requested start row is too big.

Here is an example of a request which would generate that:

http://api.geonames.org/search?featureCode=PPL&username=genadinik&country=US&style=full&adminCode1=VA&maxRows=1000&startRow=6000

How do people usually get around this?

Thanks!


texts: Rasterio resampling image to lower resolution error: 'float' object cannot be interpreted as an integer

texts:
I am trying to downsample a raster image to a lower resolution using rasterio. I am running into an error about 'float' object cannot be interpreted as an integer and I can't seem to figure out the cause of the issue or how to resolve it. I was hoping someone could suggest how to fix this.

So here is the code that I am using

with rasterio.open('image.tif') as dataset:
    data = dataset.read(
        out_shape=(dataset.height / 3.75, dataset.width / 3.75, dataset.count),
        resampling=Resampling.cubic
    )

So this error came up then. TypeError: 'float' object cannot be interpreted as an integer.

I also tried to use an integer value as well.

with rasterio.open('../openaerialmaps_imagery/ghana-alog/5bb9323e9ed15b0006d24f34-accra - ghana 2018-10-06.tif') as dataset:
    data = dataset.read(
        out_shape=(dataset.height / 4, dataset.width / 4, dataset.count),
        resampling=Resampling.cubic
    )

So I tried using a factor of 3.75(float) and a factor of 4(integer), but got the same error.

The dataset's datatype is identified as ('uint8', 'uint8', 'uint8'). But I am not sure where to go on this issue?


texts: How to Intersect a layer with custom CRS with NED data?

texts:
I am trying to create a map with a layer of county tax assessor property parcels plus NED elevation data.

The projection of the tax assessor layer is a weirdo -- looks like user-created by the county. I got the NED layer to change to that weirdo projection, or so I thought.

But when I go to intersect the NED layer with the tax assessor layer, I always get the "CRS of the layers did not match" warning, and the resultant new layer is empty.

I've tried using GRASS to re-project the NED layer, but cannot figure out how to do it that way either.

What should I do? I just want to use the parcels layer to intersect (like clip) the NED layer.

Update:

The tax assessor parcels layer gives this as its projection:

  +proj=lcc +lat_1=37.06666666666667 +lat_2=38.43333333333333 +lat_0=36.5 +lon_0=-120.5 +x_0=2000000 +y_0=500000.0000000001 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs

In the "Project Properties" window, it is listed under the heading "User Defined Coordinate Systems". What is it? And how can I get it to jive with the USGS NED data (which is, I think, unprojected and based on datum NAD83)?

Is "projection" in this QGIS usage the same as CRS?


texts:
I am having a problem using the QGIS time manager plugin when adding a raster layer. It tells me that the format of the timestamp field (I converted it to a text column) that I am referencing did not succeed, but it looks like the data is in the correct format. My column is a text type in the table, as mentioned. Is there something that you can see that I may be doing wrong (print screens below)? I was also trying to follow advice for adding the timestamp column with formatting that I found here, but I still receive an error:

https://github.com/anitagraser/TimeManager/issues/159

the last character in the next screen is my cursor, not a pipe :)

OpenLayers preview of the timestamp column (text type) from GeoServer:

QGIS 2.14.20-Essen

Ubuntu 16.04


texts: How to count the number of points in a polygon by a grouping variable in R?

texts:
I am interested in calculating the proportion of GPS points that fall within protected areas. I can do this with my full data set with the code I show below.

However, I would like to group this data by the id column and get a measure of overlap for A and B separately.

Thanks

# load packages
library(wdpar)
library(dplyr)
library(sf)

# download protected area data
# (excluding areas represented as point localities)
sa_raw_pa_data <- wdpa_fetch(c("South Africa"))
sa_raw_pa_data

# reproject data
sa_pa_data <- st_transform(sa_raw_pa_data, 4326)

# plot it
# plot(sa_pa_data)

# create some random data
set.seed(123)
d = data.frame(long=runif(100,16,38), lat=runif(100,-47,-22))
d$id <- c(rep("A", 50), rep("B", 50))

# measure the overlap
pts <- st_as_sf (d, coords = c ("long", "lat"), crs = 4326)
res <- st_intersects (pts, sa_pa_data)

length (unlist (res)) / nrow (pts) # fraction of intersecting points


texts: How to increase mask buffer and then offset it using r.mapcalc in GRASS?

texts:
I'm in the process of learning on how to remove clouds and cloud-shadows from Landsat images using GRASS. So far I was able to remove a lot of clouds from those images using i.landsat.acca, but it seems like they still require a lot of refinements:

I followed the GRASS manual for i.landsat.acca to obtain the image above: http://grass.osgeo.org/manuals/html70_user/i.landsat.acca.html.

Basically, the process is to run different Landsat band images for the same scene with i.landsat.toar, and then processed with i.landsat.acca to get the cloud areas. Then use the cloud area as a mask using r.mapcalc:

r.mapcalc "MASK = if(isnull(clouds))"

Problem#1: the cloud edges is still visible. I need to increase the cloud-mask buffer to catch those renegade cloud edges.

Problem#2: the cloud-shadows is still visible. If I include the -s flag to include a category for cloud-shadow when running i.landsat.acca, a large amounts of shadow from the terrain will consequently be removed as well... So to catch the cloud-shadows, I believe I need to make a new mask for them - by offsetting the current cloud mask until they cover the cloud-shadow areas.

To achieve both objectives, I believe I need to modify and play with the r.mapcalc command, but I can't figure out how to write it (due to my bad Math and my-even-worse programming skills)..

So, my ultimate question:

Question#1: How do I edit r.map.calc to increase the mask area (let's say by a buffer of 10 units)?
Question#2: How do we offset a mask in r.map.calc (in this case a little to the south-west direction)?

If I know the answer to both of these: I can combine both cloud-mask and cloud-shadow-mask, to get a relatively cloud-free image. Then it's a matter of overlaying it to images of the same scene of different dates.

I appreciate any help that I can get. Many thanks!


texts: Add new points into a shapefile iteratively with PyQGIS

texts:
I would like to create a filtered layer with the same attributes than an other layer, based on a condition. Moreover, I would like to fill in the shapefile iteratively when my condition is verified.

I tried to do so but my new shapefile contains only one point, when it should contain more. How I can fix it to get the whole lines into one new shp layer?

This part of script below is contained into :

For value1 in DatesHours:
 [...]

 For value2 in RefValues :
  [...]

  writer1= QgsVectorFileWriter(dir + "/Selection1.shp", "CP1250", provider.fields(), provider.geometryType(), provider.crs(), "ESRI Shapefile")
  writer2 = QgsVectorFileWriter(dir +"/Selection2b.shp", "CP1250", provider.fields(), provider.geometryType(), provider.crs(), "ESRI Shapefile")
  outelem=QgsFeature()

  for feature in macouche.getFeatures():

   true=0
   attri1=feature.attributes()

   if attri1[0] == value2 :
    ecart= datetime.datetime.strptime(value1,'%d/%m/%Y %H:%M') - datetime.datetime.strptime(feature['DATE3'],'%d/%m/%Y %H:%M')

    if ecart < zero:
     ecart= datetime.datetime.strptime(feature['DATE3'],'%d/%m/%Y %H:%M') - datetime.datetime.strptime(value1,'%d/%m/%Y %H:%M')
     if ecart_minN > ecart :
      ecart_minN=ecart
      DHp=feature['DATE3']
      true = 1
    else:
     if ecart_minP > ecart :
      ecart_minP=ecart
      DHn=feature['DATE3']
      true = 1

   if true == 1 :
    outelem.setGeometry(feature.geometry())
    outelem.setAttributes(feature.attributes())
    writer1.addFeature(write(outelem,feature))


texts: Python GUI for QGIS plugin ; how to reference form controls in code?

texts:
I've been going around in circles with this for a few hours so hopefully someone can help.  I am familiar with VB, C# but not Python - and struggling to do something fairly simple.

I have created a form using QT Designer.  On this form are a QlistView widget and a "Browse..." button.  I want to be able to select files by clicking on the button, and list these in the QlistView.  Should be easy enough.

In languages I am familiar with I would reference the QlistView object using syntax like [form name].[control name] but cannot do that here, apparently.

I have got self.dlg set to the form itself.

The list widget is a QlistWidget, with the name listWidget.

I am using the following code to handle the button click:

def handleBrowseButton(self):
pth =os.path.normpath('E:')
filter = "TXT (*.txt)"
file_name = QtGui.QFileDialog()
file_name.setFileMode(3) # 3=existing files
names = file_name.getOpenFileNames(None, "Open files", pth, filter)
if not names:
  print "CANCELLED"
  return

l = self.dlg.findChild(listWidget)

...but getting a bit stuck at the last bit.  Do I have to use findChild / findChildren / something else??

The above code does run, so I haven't got a problem getting to that point.  The issue is just how to set a reference to the listWidget object on the form in order to refer to it from that point on.

My intention is to set 'l' to refer to the listWidget and then do something like

l.list.additem(filename)


texts:
I have been trying to make JP2 raster files with FWTools with limited success.  FWTools will not allow for the creation of a JP2 file at a specific file size.  However, I am able to make a JP2 with the same data within ArcMap.  Is there a file size limit with FWTools?  Does anyone know of any work arounds?  I am looking to make larger sized GeoTiffs (100+ MB) into JP2 files.

Thanks...


texts: Writing function for User Defined Contrast Enhancement?

texts:
My problem is to implement a custom contrast enhancement for the raster images. Its not a linear enhancement, it is based on a logarithmic function.

I thought I can update the renderer block like we can do for the data provider. But it seems the writeBlock is available only for the data provider.

layer = iface.activeLayer()
rdp = layer.dataProvider()
block = rdb.block(1,extent, cols, rows)

This block gives the actual values from the raster. Not the ones from the renderer.

How do I write a function for custom Enhancement?


texts:
I have a GP tool that creates a PDF report from a map click. I'm testing running multiple reports simultaneously from different browsers. If I run two reports, one will return a PDF to the popup, and the other displays the PDF for a brief instant and then displays "Unable to generate report." If I click on the GP tool again, the report will popup, because I have set the script to check to see if a report has already been generated for that location. In the console log, I see "Timeout Exceeded." This happens after about a minute, even thought the report has obviously been generated. Do I need to set a higher timeout limit?

Here is my code snippet. The clearInterval function is for animated dots while the report is being generated, not related to the timeout.

function gpreportResultAvailable(results, messages){
                 domAttr.set(dom.byId("reportLink"),"innerHTML", "Forecast Report");

                 window.map.infoWindow.clearFeatures();

                 var contentResult;

                 results = true;

                 if(results){
                       contentResult ='';
                 }
                 else{
                       contentResult = "Unable To Generate Report";
                 }
                 clearInterval(dots);
                 window.map.infoWindow.resize(256,256);
                 window.map.infoWindow.setContent(contentResult);
                 window.map.infoWindow.setTitle(displayCountry + " Forecast Analysis Report");

                };


texts: Programmatically performing viewshed on points from a shapefile?

texts:
I want to perform a viewshed analysis on multiple points in a shapefile. I'm working with a solution as seen in this post by Jakob

Viewshed / Line of Sight analysis- as a batch process?

I modified his code such that I am able to pass the extents into r.viewshed. this is given by the variable "extcoordStr" in the below code

import processing
from PyQt4.QtCore import QFileInfo

rasterLayer = qgis.utils.iface.activeLayer()
      pointLayer=QgsMapLayerRegistry.instance().mapLayersByName('Tasmania_select_two')[0]

i=0
for ft in pointLayer.getFeatures():
   point = ft.geometry().asPoint()
   coordStr = '%d,%d' % (point.x(),point.y())
   coordArray = coordStr.split(',')

   #extent generation, +- 50,000Metres from point origin(coordStr)
   xMin = int(coordArray[0]) - 50000
   xMax = int(coordArray[0]) + 50000
   yMin = int(coordArray[1]) - 50000
   yMax = int(coordArray[1]) + 50000

   #Type juggling to get xmin,xmax,ymin,ymax
   extcoordStr = str(xMin) + ',' + str(xMax) + ',' + str(yMin) + ',' + str(yMax)

   outputViewshed = 'C:\Users\nightswitch57\Desktop\Script\viewshed_%i.tif' % i

   #running viewshed with observer and target heights set to 20Metres at max distance of 50Km
   outputs_0=processing.runalg("grass7:r.viewshed", rasterLayer,coordStr,20,20,50000,False, extcoordStr,0,outputViewshed)
   i = i + 1

problem
when I click run, my computer freezes up a bit task manager shows GRASS is been used but at the end there is nothing in the output folder.
in Jakobs answer he uses, I feel this might be the reason why I do not have any output in the folder I specified

outputViewshed = 'd:/temp/viewshed_%i.tif' % i

while I use

outputViewshed = 'C:\Users\nightswitch57\Desktop\Script\viewshed_%i.tif' % i

Also in addition I would much rather have the viewshed named after the point features "Site ID" which is in a column within the point shapefile (Tasmania_select_two) as currently it is supposed to name each viewshed as the number of the point feature.
to achieve this i'm guessing I would need to append the point features Site ID to the outputViewshed variable somehow but I don't know how to.


texts: Circle area calculation in QGIS is inaccurate

texts:
I created a circle with a radius of 3 km, calculated the area using field calculator and got a different result from real area.

The QGIS result is 28273078.8 while manually calculated area is 28274333.9 which is also the same number I got when I calculated the area in ArcMap.

"On the fly" projection is disabled and project CRS is set to utm32

Any ideas what is causing this, how to get more accurate area calculations?


texts: Total parkland acreage within each of Chicago's community boundaries?

texts:
I'm new to QGIS.

I'm trying to solve the following problem, using the tutorial here:
http://www.qgistutorials.com/en/docs/performing_spatial_joins.html
but haven't had any luck.

I have a QGIS file with two forms of spatial area data: Chicago's 77 community area boundaries on one layer, and the city's park system on another layer. The acreage of each park on the parks layer is in the ACRES column of the attribute table. What I would like to do, is find out what the total parkland acreage is within each of the 77 community areas. So if a community boundary line divides a park in half, it would obviously need to have half of its area contribute to the figure for community A and half of its area contribute for the parkland acreage in community B.

I have tried using Vector > Data Management Tools > Join attributes by location, with various options ticked and unticked, but every time it produces a new layer with a useless, blank attribute table. Then I thought perhaps this was the wrong way to go about it, so I did some more Googling and tried Vector > Geoprocessing Tools > Intersection, but this also produced a blank layer.

For your reference, I've zipped my folder containing the working file (10-mapping.qgs) and the layer files in a subfolder. Also included a screenshot of the two layers, Parks_Aug2012 and Community areas, turned on.
Google drive link: https://drive.google.com/file/d/0B_uEI02m4m1OLUdnTTh3blR2OWM

I'm using QGIS 2.18.3


texts: Getting coordinates of pixel from georeferenced drone photo?

texts:
I have a problem with geo-referenced photos from a DJI drone. I have the coordinates of the center of the photo and the metadata related to the photo and I would like to know how to go from a pixel to its relative coordinates.

Example metadata [Exif]:

latitude: 45.6426621388889
longitude: -94.5975223611111
AbsoluteAltitude: 332.743835
FlightPitchDegree: 0.3
FlightRollDegree: 3
FlightYawDegree: 89
GimbalPitchDegree: -90
GimbalRollDegree: 0
GimbalYawDegree: 7.8
RelativeAltitude: 40.700001
TlinearGain: 0
altitude: 332.7
altitude_ref: 0
aperture: 1.25
f_number: 1.25
focal_length: 13
image_make: "DJI"
image_model: "FLIR"

Having these and knowing that you have a box positioned at the pixels: 113 (xmin), 360 (ymin), 139 (xmax), 406 (ymax) how can I get its coordinates?

From a software I was able to extrapolate a solution that brings me back:

-94.6002156260815 (minimum longitude)
45.6413373184107 (minimum latitude)
-94.5998767997666 (maximum longitude)
45.6417408827611 (maximum latitude)

The problem is that I cannot understand whether it is right or whether it is calculated.


texts: Seeking podcasts related to geospatial technologies?

texts:
I would like to make a list of podcasts related to geospatial technologies.

Do you listen to any that you find interesting and informative?


texts: Problem projecting lat/long into the CRS of my raster in R

texts:
I am working with a raster in the .e00 format, located here:

ftp://ftp.epa.gov/castnet/tdep/grids/n_dw/n_dw-2000.zip

I unzip, and then load this raster in R:

require(raster)
require(rgdal)
require(sp)
r <- raster('/path/to/n_dw-2000.e00')

Which works fine:

> r.d.2000
class       : RasterLayer
dimensions  : 775, 1440, 1116000  (nrow, ncol, ncell)
resolution  : 4134.383, 4134.383  (x, y)
extent      : -2950369, 3003143, 115686.8, 3319834  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_1=0 +lat_2=29.5 +lat_0=45.5 +lon_0=0 +x_0=0 +y_0=-96 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
data source : /Users/colin/Documents/Talbot Lab/FIA_2016_manuscript_EM-ndep/extracting_e00.files_dec2016/dry_dep/n_dw-2000.e00
names       : n_dw.2000
values      : 0.4092157, 72.94985  (min, max)

I then load some lat/long values I wish to extract the values from the raster for.

Here are 6 sites, with columns corresponding to longitude and latitude:

points <- structure(c(-86.030373, -85.990492, -85.847991, -87.954285, -87.674475,
-85.138866, 46.414597, 46.137664, 42.258794, 46.567187, 45.730901,
44.383864), .Dim = c(6L, 2L))
> points
          [,1]     [,2]
[1,] -86.03037 46.41460
[2,] -85.99049 46.13766
[3,] -85.84799 42.25879
[4,] -87.95428 46.56719
[5,] -87.67448 45.73090
[6,] -85.13887 44.38386

I need to project these lat/long values into the coordinate reference system (CRS) of my raster. To do this I first convert the lat/long to a SpatialPoints data frame:

points.2 <- SpatialPoints(points, proj4string = CRS("+init=epsg:4326"))

I then re-project these points into the CRS of the raster using spTransform

points.3 <- spTransform(points.2, projection(r))

Finally, I would use the extract function to pull values from the raster for all of my points using:

out <- extract(r, points.3)

However this only returns NA's, to my dismay:

> out
     [,1]
[1,]   NA
[2,]   NA
[3,]   NA
[4,]   NA
[5,]   NA
[6,]   NA

I think the problem is how I am re-projecting the points, as if I plot the points vs. the raster they seem to be in a different orientation.

Furthermore, if I go to plot the points on otp of the raster I do not have any luck. I am unsure what specifically is going wrong however.


texts: How to add a large inset svg or pdf in QGIS 2.3 print composer?

texts:
I have a 48x36 wall map that will include 5 large inset maps, the largest of which is almost 11x17.  Due to a bug with html annotations, I need to export each inset individually then include them as images on the final large wall map.

I'd like to include them as vector art (svg or pdf).  svg has known clipping problems and does not look good in this case - and I want to see the grid coordinate labels so using a layer level masking scheme won't work.  So, I export a pdf which looks perfect in an external pdf viewer.  But, trying to include that pdf as an image in the composer just shows a large red X.  The implication is that pdf is not a valid format for including images in the composer - is that true?

As a third option, I try to export the inset map as a raster image file (png or such) but it won't do it because the resolution and the image size are too big and it gives a memory allocation error.  Reasonable but too bad.

Any other ideas?  The next try will be just writing an image of the html annotations and nothing else, and displaying them on top of all other layers in the composer.  Surely there's a better way...

Granted, if the html annotation bug (does not respect the 'Lock layers' setting in composer) were fixed, this would all be a moot point, but there's no response on that one currently, I think they have bigger fish to fry at the moment.

Thanks in advance for any ideas!

UPDATE: workaround is working so far but is lengthy:
1) export an svg of just the html annotations (all other map layers off) in just the inset area
2) edit that svg to set fill-opacity to 0 for the relevant objects (easier said than done - largely trial and error) so that the svg can be used as an overlay
3) add the edited svg as an image to the final wall-map composer

So far, the above is working for one inset.  Next steps to attempt:

repeat this for each of the 5 insets, and also for the main map since it has its own annotations that you don't want to show in any of the insets; then make sure to turn off all annotation layers for the final composer rendering.

It would be nice if the svg would export with the appropriate transparency so it can be used as an overlay in this manner, but, I'm not even sure how to word that as a feature request.  Can anyone help clarify?

UPDATE 2: (more notes on the workaround)
when saving the svg from each inset-only composer:

check 'background' for the map object, and use some color that will make it easier to search-and-replace lines in the svg file after export (i.e. look for lines with that fill= color and fill-opacity="1", and change to fill-opacity="0"
when adding the image in the main composer, make sure to uncheck 'Background' so you can see the effective transparency of the svg itself - otherwise you will go on a wild goose chase trying to find the right fill-opacity setting to change!
when editing the svg, search and replace by color as above, and you will also need to change fill-opacity from "1" to "0" for at least one more object, usually the fifth entry in the svg; if that's not it, use trial-and-error; can just hit enter on the filename of the image in composer to instantly reload and redisplay it
frame: do not display frame or grid (or coordinate labels) in the overlay svg; display those in the map object in the wall map composer


texts: How to select polygon close to each other?

texts:
I have grass patches and I put a buffer of 250m around them. I did a "dissolve to get larger polygons that overlap,  but the attribute table emitted a lot of information, which are important such as the hectare och area. How do I dissolve, and still keep these vital information.
Now the next step I want to accomplish is to "merge" grass patches that are less than 500m apart into a larger polygon. I'm having a very hard time doing this.

I'm using ArcGIS for Desktop 10.1.

What tools would be recommended?


texts: Neighbouring feature that share more than 2 vertices in Atlas

texts:
Is there a way to filter polygons that share more than 2 vertices with a feature?
I have a problem where I have a layer being controlled by Atlas, and I would like to show on Atlas only the coverage feature and the neighbouring features that share more than 2 vertices with it.
Up until now I can achieve what is in the first image using touches($geometry, @atlas_geometry) in the filter in rule-based symbology. But I wanted something like in the second image. Is there a way to do that?


texts:
I was studying the Trulia crime map and found that they were generating dynamic tiles with color block based on crime count. I am a newbie to mapping world and found how to generate static tiles, serve them, cache them etc. I am not able to understand how to generate dynamic tiles based on crime count.

On the data side, i have the latitude and longitude about where the crime happened. Is it possible to generate dynamic polygon tiles with this data ?


texts: Embedding a map in an "iframe", can't disable scrolling so it fills the whole viewport

texts:
I've tried three different ways to get this iframe / map to fill the viewport and the viewport alone:

html "iframe" scrolling options are not supported by html 5
css "style" declarations (see code) aren't working as I am implementing them.
javascript code (see code) isn't working as I am implementing it.

I got the bits of code from other stack exchange posts, researched them and attempted to apply them as seen below, with no results whatsoever. I'm sure it's just my lack of understanding how and where to apply these measures...

With the code, please note that I placed the javascript and the css where I thought they might go, and both are in the sample code just to illustrate what they are. I TRIED SEVERAL CONFIGURATIONS - -- - drawing a complete blank.

The page can be viewed here

Code:

<!DOCTYPE html>

    <meta charset="UTF-8">
    <meta name="embed-container" content="width=device-width, initial-scale=1">
    <title>Portland Cycle Safety Map</title>
    <link rel="icon" href="/icons/cog.svg">


texts: Accessing functions from another project in vb2010

texts:
I have written a project which has common functions in Visual Basic 2010 i.e. Visual Studio 2010. for e.g.

Public Function FlashPolygon(pDisplay As IScreenDisplay, pGeometry As IGeometry)
    Dim pFillSymbol As ISimpleFillSymbol
    Dim pSymbol As ISymbol
    Dim pRGBcolor As IRgbColor

    ' Set the Actions based on Error Mode
    If DEBUG_MODE = 0 Then On Error GoTo FlashPolygon_Error
    If DEBUG_MODE = 2 Then On Error Resume Next

    pFillSymbol = New SimpleFillSymbol
    pFillSymbol.Outline = Nothing

    pRGBcolor = New RgbColor
    pRGBcolor.Green = 128

    pSymbol = pFillSymbol
    pSymbol.ROP2 = esriRasterOpCode.esriROPNotXOrPen

    pDisplay.SetSymbol(pFillSymbol)
    pDisplay.DrawPolygon(pGeometry)
    Sleep(300)
    pDisplay.DrawPolygon(pGeometry)

    On Error GoTo 0
Exit Function

I want to access this function in a form in another project. I have already added reference to the commonfunctions project. But i am getting error Reference to a non-shared member requires an object reference.

Can someone help with what mistake I am doing in this code?


texts: Is it possible to draw a 3D Choropleth Map in ArcGIS?

texts:
What I have in mind is something like these two examples below:

Also, if it is not possible to create a 3D choropleth map in ArcGIS, is there some other GIS software package that can draw 3D choropleth maps other than Maptitude?


texts: How to split ESRI file geodatabase into two or three by region?

texts:
I have got a ESRI file geodatabase(10.3) with many feature-classes(Polygon, Polyline & point). I don't have enterprise ESRI license. I have to split the entire data into two unique parts to do simultaneous edit/update by two employees.

How can I split the file geodatabse(all feature-classes) into two portions, each containing all the feature classes as shown in the image?


texts: How to bring my PostGIS map from GeoServer to Leaflet based web application?

texts:
I have created a group of shapefiles using QGIS and I have stored their corresponding information in PostGIS database by importing the shapefile into pgAdmin3 using shapefile and DBF loader. Then I imported the whole PostGIS database to GeoServer. It works fine using OpenLayers link in layer preview option of GeoServer. The attribute values are showing up in tabular form when I click on each object.

What I need now is to create a web based application using Leaflet or OpenLayers where the attribute value is showing as pop-ups and I want the facility to add more shapefiles/PostGIS table dynamically.


texts:
I have a Geoserver (version 2.5.1) with layers that are assigned to specific users and groups which has been working very well for quite a while. Operationally very little has changed. New users & groups have been added and new layers published and assigned to the groups. Layers for each group are published to a specific group are published to a common workspace. some of the layers use a common style which is published to a shared workspace. For example, a shapefile in Workspace (X) and another published to Workspace (Y) both use a style from Workspace (Z).

This has been working well for about 2 years. Suddenly some users are not able to view any layers. When I checked the GetCapabilites XML it is being truncated at about row 6440. The error points to some error with a specific layer.

<Layer queryable="1" opaque="0">
    <Name>kuilsriver:Kalkfontein BW</Name>
    <Title>Kalkfontein Block Watch</Title>
    <Abstract/>
    <KeywordList>
      <Keyword>features</Keyword>
      <Keyword>Kalkfontein</Keyword>
    </KeywordList>
    <SRS>EPSG:4326</SRS>
    <!--<?xml version="1.0" encoding="UTF-8" standalone="no"?><!DOCTYPE
    ServiceExceptionReport SYSTEM

"http://live.webserver.com:80/geoserver/schemas/wms/1.1.1/WMS_exception_1_1_1.dtd">
      javax.xml.transform.TransformerException: java.lang.NullPointerException
java.lang.NullPointerException
null

If I delete the layer the error is repeated for the next alphabetical layer, and so on.

There have been no changes to the server except for adding users, groups, roles and layers so why this error should suddenly occur is unexplained.

I initially thought that the error may be caused something in a recently published layer but after removing all the layers published in the last week (with the corresponding users, groups and roles) the error persists.

I have checked other posts about errors in the GetCapabilites XML but resolution suggestions have not solved my problem.


texts: Open library for generating Styled Layer Descriptor (SLD) files?

texts:
I'm developing an app that can create/edit shapefiles and visual representation of vector layers and to export them to some of most common geoportals(f.ex. GeoServer).

I have dealt with .shp, so now i need to deal with styles' code.

I already looked into C++ libs for XML creation, but I think there might be a specific lib for just exactly SHP+SLD exportation.

Is there a library (preferably C++) for SLD export of vector layer styles?


texts: Is there an automated way to import a GPX file and use the filename as the layer name in QGIS?

texts:
I'm using QGIS, importing a lot of GPX files, always selecting the tracks layer for each one and clicking OK. The layers end up being named tracks by default.
Is there a way to import GPX files (individual or bulk) and automatically use the filename as the layer name instead of having to manually rename the layer each time I import a file?


texts: Plotting only raster cells that fall within boundaries of another raster in R

texts:
I've got two rasters, soil and slope. I want to plot the slope raster but only in areas where there is soil (i.e. soil is > 0). How can I do this?

This is the soil raster:

This is the slope raster:

EDIT 1: more information about the rasters.

soil

 > extent(soil)
 class      : Extent
 xmin       : -180
 xmax       : 179.9856
 ymin       : -60
 ymax       : 89.994
 > res(soil)
 [1] 0.008333 0.008333
 > origin(soil)
 [1]  0.001133 -0.002400

slope

 > extent(slope)
class      : Extent
xmin       : -180
xmax       : 180.0023
ymin       : -56.00016
ymax       : 84.00257
> res(slope)
[1] 0.008333386 0.008333496
> origin(slope)
[1] 0.0011330 0.0009332

EDIT 2: information about dimensions

 > dim(soil)
 [1] 18000 43200     1

 > dim(slope)
 [1] 16800 43200     1


texts: How to split a feature and save back to database?

texts: How to split/cut a polygon through OpenLayers split method and same editing save to the PostGIS?


texts: Loading older google satellite image with OpenLayers plugin in QGIS?

texts:
I am using the OpenLayers plugin in QGIS. I want to use the google satellite layer, but there is a cloud on the specific image.

Is it possible to load one of the older google satellite images?

Like loading the older satellite images on google earth.


texts: Changing position of polygons in ArcGIS Desktop?

texts:
I am using ArcGIS Desktop 10.4.1.

After converting the existing polygons to points, we moved the points to the desired point using Network Analyst.

Here, the FID value of the existing polygon and the moved point SHP file are the same.

I want to use this to move the coordinates of the existing polygon center to the coordinates of the moved point.

The reason for this is that I want to keep the rotation angle of the existing polygon anyway.

If you make the point back to polygon, the angle will be 0 degree in a lump.

The coordinates of the points are obtained with Arcgis's "Add x y coordinates", and I'm wondering how I can change the position of the polygons.

That is, if you define it briefly, I want to move to the coordinate value of the point, not the existing coordinates of the polygon.

Is there a way?


texts: What are the existing datasets of world maritime routes?

texts:
I have a list of cargo ports of all over the world and need to show some port-to-port trips on a map, and also to compute lengths of these trips.

For that, I need a dataset of actual maritime routes (even approximative) between important ports preferably free. Is there something available?


texts: Retrieving features from QgsFieldExpressionWidget using PyQGIS?

texts:
I'm working on a plugin for QGIS 3 and I'm fighting with the QgsFieldExpressionWidget widget.

I don't know how how to evaluate the expression added in the box and get the corresponding features as a python list.

Within the plugin I would like to write a function to call that is able to know if in the QgsFieldExpressionWidget is a simple field (without expression) or a expression that has to be evaluated.

No problem for the first part, but API changed for the QgsFieldExpression and I'm not able to follow the instruction of the pyqgis cookbook.

This is the function I'm trying to write:

 def getFields(lay, exp):

    # vector layer from the combobox
    vl = lay.currentLayer()
    # vector field (or expression)
    field = exp.currentText()

    data = []

    # get the data from the layer
    # no expression
    if not exp.currentField()[1]:
        for i in vl.getFeatures():
            data.append(i[field])

    # if expression is selected

    # this is the piece of code that is not working for QGIS 3

    else:
        # expression name
        fil = exp.currentField()[0]
        # expr = QgsExpression(fil)
        # for i in vl.getFeatures():
            # data.append(expr.evaluate(i, vl.pendingFields()))

    return data


texts:
In many Earth Engine functions such as Image.reduceRegions or Image.sample, there is an argument called tileScale. The API documentation generally interprets it as:

A scaling factor used to reduce aggregation tile size; using a larger
tileScale (e.g. 2 or 4) may enable computations that run out of memory
with the default.

While this seems to be well documented, I still do not understand the effect of tileScale on the final result of the computation.

What does tileScale do?
Does tileScale chunks the image before performing the computation?
Does it have an effect on the scale?

Below is a code sample to work with (code link).
// Getting the image of the region of interest
var roi = ee.Geometry.Point([1.864578244475683, 14.492292970253338]);
var image = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
                .filterDate('2019-01-01', '2019-01-31')
                .filterBounds(roi)
                .select(['B5', 'B4', 'B3'])
                .toBands()
                .rename(['B5', 'B4', 'B3']);

// Checking it out
print(image);

// Define the visualization parameters.
var vizParams = {
  bands: ['B5', 'B4', 'B3'],
  min: 0,
  max: 0.5,
  gamma: [0.95, 1.1, 1]
};

// Center the map and display the image.
Map.centerObject(image, 9);
Map.addLayer(image, vizParams, 'false color composite');

// Computing the band means at 30 meters and tileScale of 1: (Computation succeeded)
var control = image.reduceRegion({
  reducer:ee.Reducer.mean(),
  scale    :  30,
  tileScale:   1,
  maxPixels:1e13
});

print(control, 'first scenario');

// Computing the band means at 5 meters and tileScale of 1: (Computation succeeded)
var scenario2 = image.reduceRegion({
  reducer:ee.Reducer.mean(),
  scale:        5,
  tileScale:    1,
  maxPixels: 1e13
});

print(scenario2, 'second scenario');

// Computing the band means at 5 meters and tileScale of 16: (Computation Error: Computation timed out.)
var scenario3 = image.reduceRegion({
  reducer:ee.Reducer.mean(),
  scale:        5,
  tileScale:   16,
  maxPixels: 1e13
});

print(scenario3, 'third scenario');


texts: Setting up a relationship class between two feature classes

texts:
I am working on a project that involves a line feature class that represents roads and a point feature class that represents defects on the roads. I am trying to relate the point feature class to the line feature class so that when I select the road I am able to add a related point with the same name that represents a defect on that road.

I am having no trouble when I create a table and relate records to the line feature class but when I create the point feature class and try to relate it and add a new point in ArcMap I am not able to.

Has anyone experienced this before?


texts: Achieving the opposite of the shortest_line() function in QGIS?

texts:
I currently use the Geometry Generator styling on a point layer to draw the intersection between protection buffers around tree points and another polygon (say a building footprint) using the following expression:

difference(buffer($geometry,"buffer"),
           aggregate('building_footprint_layer','collect',$geometry))

The style typically calls for a line from the tree point to some edge of the polygon to make it clear which tree is connected to which polygon.

Currently I can achieve this on these irregular polygons using shortest_line(), which is useful (see blue lines in image below).

(intersection() would be another option per another question of mine)

But I would also like to draw a line from the tree point to the furthest edge of the resultant polygon so I can slap a label on it and indicate what the protection buffer radius is (outside of the intersection with the building polygon). See black lines in image above.

Keeping in mind that

The point will always be within the buffer polygon
Intersections with the building polygon can be at any angle so I can't just draw a line at say 45° and take the intersection of that.
The final angle of the longest line to the polygon doesn't matter much

Is there any way to achieve this using expressions and the Geometry Generator style? I have a feeling a custom Python function would be required; my understanding is that finding the longest linear distance from a point to a polygon boundary would require some kind of sweeping algorithm (?)


texts: Using Field Value as Search Distance Parameter in Select by Location in ArcPy?

texts:
I've been trying to use the select layer by location tool in a loop to select features that are within a specific distance to one feature class for each feature in a different feature class (each retaining wall has a search distance around it of twice its height). Unfortunately I've realized from Esri Support that it isn't possible to use a FIELD as the argument for the "Search Distance" Parameter.

Is there any other way to do what I'm trying to do here without using the Select Layer by Location tool? I'll post the relevant pieces of my code below.

Retain=arcpy.GetParameterAsText(0)
TSafetyRW=arcpy.GetParameterAsText(1)

if Retain=="":
    Retain=dataPath+"/LayerHome.gdb/Retaining_Walls"
if TSafetyRW =="":
    TSafetyRW=dataPath+"/LayerHome.gdb/Transportation_Safety_RW"

inRetain=os.path.splitext(os.path.basename(Retain))[0]
inTSafetyRW=os.path.splitext(os.path.basename(TSafetyRW))[0]

TranspoSafetyScore = "Socio_Transpo_Safety"
Wall_Height = "RW_Height"
Road_to_Wall_Dist = "Distance"
Road_Class = "Class_Type"
Twice_Height = "Threshold_Distance"

##Adding and assigning a default value to some fields
arcpy.AddField_management(inRetain, TranspoSafetyScore, "FLOAT")
arcpy.CalculateField_management(inRetain, TranspoSafetyScore, 1)

##Tuples used to assign a Road Class to a Retaining wall
Road_to_Wall = [("No Road", 1),
                ("MUP", 1),
                ("Sidewalk", 1),
                ("Collector", 2),
                ("Major Collector", 2),
                ("Railway", 3),
                ("Arterial", 4),
                ("Transitway", 4),
                ("City Freeway", 4)]

##Tuple calculation which creats a temporary layer of the TSafetyRW FC for each road class
for pairs in Road_to_Wall:
    (roadclass, score) = pairs
    arcpy.MakeFeatureLayer_management(TSafetyRW, "justselected", "\""+Road_Class+"\"='"+str(roadclass+"'"))
    arcpy.SelectLayerByLocation_management(inRetain,"WITHIN_A_DISTANCE","justselected", Twice_Height,"NEW_SELECTION")
    arcpy.CalculateField_management(inRetain, TranspoSafetyScore, score)
    arcpy.SelectLayerByAttribute_management(inRetain, "CLEAR_SELECTION")
    arcpy.Delete_management("justselected")

As I mentioned above, I know that you can't have a field as the search distance parameter for the Select by Location tool so I know that this code won't run properly.

I specified the type of variable that I wanted to use for the search distance parameter in the first paragraph

In other words, I want to use a field from the Target Feature Class as a variable, so it's constantly changing with each retaining wall as the loop is executed. In this case I want to use the Twice_Height variable as the search distance, which is the field "Threshold_Distance" which is a part of the target Feature Class "inRetain".


texts:
I receive "mapserv request processing time (msLoadMap not incl.)" in my ms_error.log, but no idea what that means, is that an error or the normal behavior of MapServer?

It appears regardless the layer is point, line or polygon:

Point example:

[Mon Oct 27 07:56:34 2014].499867 CGI Request 1 on process 7956
[Mon Oct 27 07:56:34 2014].505409 msDrawMap(): rendering using outputformat named png (AGG/PNG).
[Mon Oct 27 07:56:35 2014].480974 msDrawMap(): Layer 0 (Boreholes), 0.976s
[Mon Oct 27 07:56:35 2014].481018 msDrawMap(): Drawing Label Cache, 0.000s
[Mon Oct 27 07:56:35 2014].481032 msDrawMap() total time: 0.980s
[Mon Oct 27 07:56:35 2014].573357 msSaveImage(stdout) total time: 0.092s
[Mon Oct 27 07:56:35 2014].574595 mapserv request processing time (msLoadMap not incl.): 1.075s
[Mon Oct 27 07:56:35 2014].574638 msFreeMap(): freeing map at 0x80012a68.

Line example:

[Mon Oct 27 07:56:36 2014].69162 CGI Request 1 on process 8928
[Mon Oct 27 07:56:36 2014].84536 msDrawMap(): rendering using outputformat named png (AGG/PNG).
[Mon Oct 27 07:56:37 2014].649376 msDrawMap(): Layer 7 (Seismic_profiles), 1.565s
[Mon Oct 27 07:56:37 2014].649442 msDrawMap(): Drawing Label Cache, 0.000s
[Mon Oct 27 07:56:37 2014].649455 msDrawMap() total time: 1.579s
[Mon Oct 27 07:56:37 2014].706250 msSaveImage(stdout) total time: 0.057s
[Mon Oct 27 07:56:37 2014].707153 mapserv request processing time (msLoadMap not incl.): 1.638s
[Mon Oct 27 07:56:37 2014].707185 msFreeMap(): freeing map at 0x80012a80.

Polygon example:

[Mon Oct 27 07:56:36 2014].1015 CGI Request 1 on process 4152
[Mon Oct 27 07:56:36 2014].6791 msDrawMap(): rendering using outputformat named png (AGG/PNG).
[Mon Oct 27 07:56:36 2014].245940 msDrawMap(): Layer 8 (Forearc_basins), 0.239s
[Mon Oct 27 07:56:36 2014].245976 msDrawMap(): Drawing Label Cache, 0.000s
[Mon Oct 27 07:56:36 2014].245986 msDrawMap() total time: 0.244s
[Mon Oct 27 07:56:36 2014].297990 msSaveImage(stdout) total time: 0.052s
[Mon Oct 27 07:56:36 2014].299173 mapserv request processing time (msLoadMap not incl.): 0.298s
[Mon Oct 27 07:56:36 2014].299201 msFreeMap(): freeing map at 0x80012a70.

Does "msLoadMap not incl." represent an error or a bottleneck problem because of cygwin?

OS:
win8.1 i5 (1.6GHz,8GB RAM), mapserver under cygwin

$ mapserv -v
MapServer version 6.4.1 OUTPUT=PNG OUTPUT=JPEG SUPPORTS=PROJ SUPPORTS=AGG SUPPORTS=FREETYPE SUPPORTS=CAIRO SUPPORTS=ICONV SUPPORTS=FRIBIDI SUPPORTS=WMS_SERVER SUPPORTS=FASTCGI SUPPORTS=GEOS INPUT=JPEG INPUT=POSTGIS INPUT=SHAPEFILE

$ uname -a
CYGWIN_NT-6.3-WOW64 gery 1.7.32(0.274/5/3) 2014-08-13 23:03 i686 Cygwin


texts: what is the Parcel (APN) in this data block?

texts:
I am trying to figure out the 'pattern' to obtain APN to polygon mapping of   parcel in datasets...so my code can query by "APN" for a polygon.

APN as in: https://en.wikipedia.org/wiki/Assessors_parcel_number

  "An assessor's parcel number, or APN, is a number assigned to parcels
  of real property by the tax assessor of a particular jurisdiction for
  purposes of identification and record-keeping. The assigned number is
  unique within the particular jurisdiction, and may conform to certain
  formatting standards that convey basic identifying information such as
  the property type or location within the plot map."

I downloaded Balitmore,MD parcel:
https://data.baltimorecity.gov/Geographic/Parcels/rb22-mgti

When clicking a polygon I get:

bldg_no:4511
blocklot:8455 001
blocknum:8455
capture_me:BLK
cntgraph:1/2
edit_date:1997-03-18T08:00:00.000Z
fulladdr:4511 WINDSOR MILL ROAD
last_org:100
last_user:
parcelnum:001
pin:8455001
shape_area:1789137.20583
shape_leng:6749.7986708
st_name:WINDSOR MILL
st_type:ROAD
stdirpre:
subtype_ge:1
tag:600000000263749
unit_num:

What column is considered the APN?


texts:
I'm looking solution to read DGN file in java or c# without any success.
Gdal makes many problem and doesn't support all element eg. complex element.

Does anyone know of a library that supports the reading of all the elements DGN file?

Where I can find format dgn file?


texts: Obtaining a contour as a single polygon in ArcGIS?

texts:
I have a DEM that I want to retrieve a single contour polygon from at a given elevation. I can use the Contour tool to obtain a series of polylines at a given elevation by setting the contour interval to a point where only one contour is produced (e.g. 9999m) at the base contour interval (e.g. 105m).

What I want is a polygon representing the specified elevation, which happens to be the average water surface elevation of a river. You can see an example of my data in the image below.

Any thoughts on how I can accomplish this? The polyline contours run up to the edge of the DEM and then break off, so there is no continuous line which surrounds the whole river.

The second image is essentially what I would like to be output. You can see that I'll also have to solve the problem of islands within the river channel. I am attempting to automate/model this task so ArcObjects are fine.